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INTRODUCTION  AND  SIGNIFICANCE  OF  PROBLEM 


During  the  past  two  decades,  a  significant  amount  of  effort  has  been  directed 
toward  developing  methodologies  for  turbulence  closure  of  the  time-averaged  Navier- 
Stokes  (N-S)  equations.  Several  workshops  have  been  held  to  address  issues 
concerning  turbulence  modeling.  These  issues  consisted  of  the  use  of  simple  models 
to  predict  complex  flows  that  are  three-dimensional  in  nature,  computational  efficiency 
of  higher  order  models,  test  database  for  model  validation,  and  the  degree  of  accuracy 
of  higher  order  models  over  simple  models  to  predict  complex  flow  experienced  in 
engineering  design.  The  most  noted  consisted  of  the  AFSOR/HTTM-Stanford 
Conference  (Kline,  et  al,1982)  on  complex  turbulent  flows.  This  conference  followed 
an  earlier  conference  (1968)  that  extended  the  limited  two-dimensional, 
incompressible  work  to  include  compressibility,  curvature,  mass  transfer,  separation, 
and  secondary  motion.  More  recently,  several  workshops  have  been  held  on  the 
subject  matter  relative  to  application  to  the  National  AeroSpace  Plane  (NASP) 
program  (NASP,  1989  -  1991)  as  well  as  Cornell  University  (Lumley,  1990)  that  led  to  a 
continuation  of  the  Stanford  effort  (Bradshaw,  et  al,  1-91).  This  recent  work  considered 
a  more  complete  range  of  test  cases  with  a  focus  on  compressible  flows  that  included 
shock/boundary  layer  interactions. 

In  general,  the  techniques  that  have  been  developed  to  date  have  been  based 
on  two-dimensional,  incompressible,  fundamental  parallel  shear  flows.  These 
techniques  have  subsequently  been  extended  to  three-dimensional,  compressible, 
complex  flowfields  and  have  been  shown  to  be  inadequate  and  usually  without 
generality.  While  advances  have  been  made  in  numerical  algorithms  to  aid  the 
turbulence  modeler  and  facilitate  computer  storage,  full  Reynolds  stress  solutions  is 
still  considered  long  range.  Accordingly,  generalization  of  existing  techniques  based 
on  low-speed  turbulence  models  and  physical  evidence  appears  capable  of  meeting 
near  term  requirements  for  the  design  process  and  will  be  the  focus  of  this  review. 

Several  recent  papers  have  appeared  in  the  literature  that  have  essentially 
reviewed  the  state-of-the-art  in  turbulence  modeling  and  have  formed  the  basis  of  this 
investigation.  These  include  the  work  of  Rodi  (1-91)  and  Chen  and  Patel  (1988)  who 

introduced  a  hybrid  model  that  consists  of  the  standard  k-  e  equations  for  the  outer  part 
of  the  boundary  layer  and  a  one-equation  model  with  standard  empirical  description  of 
the  length  scale  (e)  for  the  near  wall  region. 

Compressibility  effects  were  examined  by  Rubesin  (1-89),  Viegas  and  Rubesin 
(6-91),  Cousteix,  et  al  (1-91),  Situ  and  Schetz  (6-91),  and  Dash  (6-91). 
Compressibility  corrections  appeared  to  improve  the  spreadrate  predictions  of 
unconfined  flows  but  did  not  work  as  well  as  turbulence  models  without  the  corrections 
for  boundary  layer  flows. 

A  critical  evaluation  of  two-equation  models  for  near  wall  turbulence  was 
presented  by  Speziale,  et  al  (6-90)  and  an  excellent  examination  of  several  turbulence 
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models  was  made  by  Menter  (June,  1991)  for  attached  and  separated  flows.  In  the 
latter  paper,  the  models  of  Baldwin  and  Barth  (1-91  and  8-90),  Johnson  and  King  (11- 
85),  Baldwin  and  Lomax  (1978),  and  Wjlcox  (11-88)  were  examined.  Wilcox  (1-91  and 

6-91)  presents  an  historical  review  of  the  k-co  two-  equation  model  as  well  as  progress 
in  hypersonic  turbulence  modeling.  The  later  paper  also  included  a  critique  of 
compressible  effects. 

Horstman  (6-91)  compared  several  hypersonic  flow  cases  using  a  standard  k-e 

model,  the  Rodi  (1-91)  hybrid  model,  and  the  k-e  model  with  variations  of  compressible 
effects  using  the  work  of  Rubesin  (6-90)  and  Verong  and  Coakley  (1-87)  for  correction 
to  the  turbulent  length  scale.  To  examine  the  effects  of  the  individual  modifications  in 
the  compressible  model,  Horstman  provided  solutions  with  the  Rubesin  correction 

alone  and  with  the  length  scale  correction  alone.  When  compared  to  the  standard  k-e 
compressible  model,  the  length  scale  correction  had  a  significant  impact  on  heat 
transfer  in  the  reattachment  region  as  well  as  controlling  the  separation  length. 


Turbulence  models  that  have  received  wide  spread  attention  are  the  k-e  models 
which  were  developed  for  incompressible  flows.  The  apparent  limitation  in  the  k-e 
models  lies  with  the  e  transport  equations  intrinsic  behavior  approaching  solid 
boundaries  (viscous,  near  wall  region).  Accordingly,  bridging  techniques  (wall 
functions)  have  been  adapted  to  accommodate  the  viscous  thin-layer  region  where 
steep  gradients  exist.  However,  the  use  of  the  law  of  the  wall  and  the  assumption  that 
turbulent  quantities  can  be  expressed  in  terms  of  the  wall  shear  velocity  fail  under 
severe  pressure  gradient  flows  since  the  near  wall  region  is  not  controlled  by  wall 
shear  stress. 

An  examination  of  two-equation  models  with  low  Reynolds  number  techniques 
was  made  by  Patel  et  al,  (1985)  for  incompressible  flows  which  indicated  that 
improvements  were  required  for  damping  in  order  to  develop  a  more  generalized 
turbulence  model.  Shih  (8-90)  reviewed  this  work  together  with  direct  simulation 
(DNS)  of  the  Reynolds  stresses  (Mansour  et  al,  1988)  and  proposed  a  modification  of 
the  pressure  transport  term  of  the  k-equation  as  well  as  a  length  scale  correction  to 
account  for  the  asymptotic  behavior  of  the  dissipation  rate  to  be  e  -  0  (y2)  near  the  wall. 

Table  I  provides  a  hierarchy  of  turbulence  models  which  has  expanded  upon 
the  summaries  of  Dash  (6-91)  as  well  as  Rubesin  (1-89).  Dash  noted  that  with  current- 
generation  main-frame  computers  (CRAY-2,  for  example),  a  grid-resolved,  steady  3D 
flowfield  solution  using  a  two-equation  turbulence  model  in  an  efficient  PNS  solver 
can  require  over  twenty-four  hours  of  CPU  time,  depending  on  grid  requirements  and 
complexity  of  thermochemistry.  Accordingly  a  practical  viewpoint  is  required  to  meet 
near-term  work  required  for  use  in  advanced  CFD  codes.  Turbulence  modeling  at  the 
two-equation  level  meets  this  requirement. 
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TABLE  1  HIERARCHY  OF  TURBULENCE  MODELS  (CONTINUED) 
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As  noted  in  Table  I,  direct  numerical  simulation  (solves  the  full/non-averaged 
Navier-Stokes  equation)  and  large  eddy  simulation  represent  the  most  complete 
physics  for  resolving  flow  characteristics;  however,  these  methodologies  are  restricted 
to  simulation  of  simplistic  shear  flows  at  low  Reynolds  numbers.  Second-order  closure 
have  not  demonstrated  invariance  with  fixed  coefficients,  and  their  inclusion  into 
advanced  computer  codes  would  not  permit  solutions  of  realistic  flow  problems  on 
current  generation  computers. 

Reynolds  stress  and  generalized  extensions  of  two-equation  models  (such  as 
algebraic  Reynolds  stress  rnodel^  have  eliminated  the  Boussinesq/isotropic 
assumptions  but  still  require  several  partial  differential  equations  for  the  stress  tensor, 
or  require  multiple  scales,  or  modifications  to  closure  constants  to  allow  for  more 
physics.  The  degree  of  complexity  of  these  models  together  with  computational  time 
have  not  shown  to  be  significantly  better  than  the  simpler  two-equation  models. 

The  remainder  of  Table  I  shows  the  standard  two-equation  models  currently 
used  in  the  scientific  community,  variations  of  these  models  that  have  been  adapted  to 
handle  near-wall  conditions,  compressibility,  and  curvature  as  well  as  several  recent 
lower  order  models  that  have  shown  successes  for  a  variety  of  complex  flows.  As 
previously  noted,  the  most  popular  two-equation  model  is  the  k-e  model  where  ad-hoc 
extensions  have  been  made  to  include  compressibility  and  curvature.  In  order  to 
eliminate  numerical  difficulties  of  integration  (of  dissipation)  to  the  wall,  bridging 
techniques  (van  Driest),  wall-functions  (low-Reynolds  number  techniques)  and  one- 
equation  models  (hybrid  models)  have  been  adopted  to  accommodate  length  scales 
for  the  near-wall  region.  In  the  latter  models,  the  use  of  kinetic  energy  has  been  used 
for  the  velocity  scale. 

A  one  and  one-half  equation  model,  so  named  in  as  much  as  it  retains 
characteristics  of  the  field  equations  (k-e  ),  has  recently  been  developed  (Baldwin  and 
Barth,  1-91)  that  couples  the  k-e  equations  through  a  transformed  variable  (Rj  = 

k2/ue)  that  eliminates  need  for  a  length  scale  (but  introduces  damping  similar  to  the 
hybrid  models  above).  Finally,  one-equation  models  that  account  for  lag  effects  and 
use  the  kinetic  energy  field  equation  for  the  velocity  scale,  variations  that  include  two- 
time-scales,  a  one-half  model  that  uses  an  ordinary  differential  equation  to  introduce 
lag  between  turbulence  and  mean  motion,  and  the  popular  algebraic  eddy  viscosity 
models  complete  the  table. 

While  many  successes  have  been  shown  for  unit  problems,  i.e.,  individual  effects 
of  pressure  gradient,  mass  transfer,  roughness,  non-isoenergetic  flow  etc.,  within  the 
incompressible  or  compressible  flow  state,  coupling  of  these  effects  has  not  been 
demonstrated.  A  turbulence  model  utilized  for  a  generalized  flow  problem  must  be 
able  to  produce  the  "simpler"  unit  level  problem  as  well  as  the  coupling  effects  in  order 
to  fully  embody  the  physics  in  the  model  and  provide  agreement  with  data  that  would 
not  be  fortuitous  or  due  to  error  cancellation.  In  attempting  to  meet  this  objective, 
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consideration  to  Reynolds  stress  models  featuring  algebraic,  one-  equation,  and  two- 
equation  models  will  be  made  relative  to  velocity  and  length  scale  characteristics. 

1 . 1  Program  Objectives 

The  objectives  of  the  proposed  program  is  to  investigate  the  embodied  physics 
of  several  popular  turbulence  closure  models  in  order  to  determine  the  feasibility  of 
developing  universality  capable  of  providing  near  term  solutions  to  engineering 
design  problems.  The  model  should  include  compressibility,  wall  conditions  such  as 
roughness,  mass  transfer,  and  pressure  gradient/curvature  and  be  applicable  to  both 
the  boundary  layer  and  free  shear-layer  flows.  While  the  focus  is  on  parallel  shear- 
layer  flows,  a  secondary  objective  was  to  examine  flows  experiencing  separation  and 
interactions.  The  module  selected  with  appropriate  modifications  should  be  adaptable 
to  state  of  the  art  multizonal  Navier-Stokes  solvers  used  in  the  industrial  community.  A 
roadmap  to  achieve  these  objectives  which  can  be  accomplished  in  a  Phase  II 
program  will  be  investigated. 

1.2  Report  Organization 

The  report  has  been  organized  into  six  sections  and  appendices.  This 
introduction  (Section  1)  provides  an  introduction  of  the  turbulence  modeling 
methodologies  and  issues,  and  is  followed  by  Section  2  that  examines  several 
turbulence  models  with  the  pertinent  physics  that  delineate  reasons  why  the  models 
work  for  various  flow  conditions.  Section  3  provides  a  recommendation  for  a 
turbulence  closure  model  that  can  be  used  in  a  multizonal  PNS  solver  (Section  4)  and 
encompass  sufficient  generality  to  handle  a  large  body  of  geometric/flow  engineering 
problems.  A  road  map  is  provided  in  Section  5  that  describes  how  the  recommended 
turbulence  model  can  be  generalized  through  a  building  block  approach  to  handle 
unit  level  as  well  as  coupling  of  unit  level  problems  subject  to  compressible  flow 
conditions.  A  compressible  flow  database  will  be  outlined  for  this  approach.  Finally, 
Section  6  assesses  the  results  of  the  Phase  I  investigation  as  well  as  outlines  a 
technical  approach  for  a  Phase  II  program  to  obtain  a  generalized  turbulence  model 
capable  of  solving  a  large  class  of  practical  engineering  problems. 

Appendices  are  also  provided  that  list  details  of  the  turbulence  models 
examined  during  this  program. 
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2.0  TURBULENCE  MODEL  REVIEW 


This  review  contains  subject  matter  that  can  be  found  in  a  number  of  texts  (or 
papers)  but  is  repeated  for  completeness  in  order  to  provide  background  concerning 
the  turbulence  scales  of  length  and  velocity  and  why  variations  in  the  choice  of  these 
scales  can  solve  a  number  of  turbulent  flow  states. 


The  review  focussed  on  models  that  employ  the  Boussinesq  criteria  that  the 
stress-strain  law  for  time-averaged  turbulent  flow  is  analogous  to  viscous  fluids.  The 
eddy  viscosity  function  can  be  expressed  in  terms  of  length  (i)  and  velocity  (u)  scales 
such  as 

Pt=  iu 

The  length  and  velocity  scales  can  be  obtained  via  algebraic  functions 
involving  properties  of  the  mean  flow  and  geometry  or  from  differential  equations  for 
one  or  more  properties  of  the  turbulent  motion.  When  i  and  u  are  determined  in  terms 
of  an  algebraic  form  of  the  mean  motion,  the  models  are  referred  to  as  zero-equation. 

When  i  is  determined  algebraically  and  u  is  found  from  a  solution  to  the  turbulent 
kinetic  energy  transport  equation,  the  model  is  referred  to  as  a  one-equation  model. 
When  both  i  and  u  are  determined  from  transport  equations,  the  model  is  considered  a 
two-equation  model 

The  Boussinesq  criteria  allows  for 

T  =  p  u'  v'  =  fit  du/dy 

where  (it  is  a  property  of  the  local  state  of  turbulence  (becomes  effective  only  when 

there  is  motion)  and  is  not  a  property  of  the  fluid.  Prandtl  proposed  a  formulation  for 
the  eddy  viscosity  (mixing  length  hypothesis),  namely 


|3y| 

where  i  the  mixing-length,  is  algebraically  prescribed  and  is  characteristic  of  the  mean 
free  path  of  molecules  (kinetic  theory  of  gases).  The  eddy  viscosity  implies  a  gradient- 
type  diffusion. 

Prandtl  later  proposed  a  model  that  employed  the  turbulence  kinetic  energy  as 
the  characteristic  velocity  scale,  such  that 

Pt  =  p  k1/2 1 
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where  i  is  again  prescribed  algebraically  and  k  requires  a  solution  of  a  differential 
equation  of  turbulent  energy  transport.  Kolmogorov  proposed  two  independent 
properties  to  represent  the  eddy  viscosity  (see  Wilcox  description,  f-91)  that  consisted 
of  turbulence  energy  k  and  a  characteristic  frequency  of  the  motion.  For  this  situation, 
the  eddy  viscosity  is  written  as 

=  p  k/f 


where  comparison  to  the  Prandtl  model  allows  for  i  =  k1/2/f  and  both  k  and  f  (which  is  co 
in  the  k-  to  models)  are  determined  by  differential  equations  of  transport. 


Several  modifications  have  been  made  to  these  basic  models.  One  of  the  more 
notable  modifications  is  due  to  van  Driest  who  introduced  a  damping  factor,  such  that 


Heff  =  |lL  +  Ht  =  m.  +  pt2m 


du 

By 


where 


im  =  ky[l-exp(-y+/A+)] 

for  u+  =  uzy/v  and  A+  =  26.  An  effective  viscosity  is  defined  that  allows  for  both 
laminar  and  turbulent  contributions  where  the  latter  includes  a  damping  function 
whose  effect  is  experienced  in  the  near  wall  region  (y+  <  A+). 

The  above  represent  a  basic  formulation  from  which  eddy  viscosity  models  are 
developed.  What  follows  will  be  a  discussion  on  variations  of  length  and  velocity 
scales  that  have  been  recently  reinvestigated  and  applied  to  a  variety  of  complex 
flows.  The  Cebeci-Smith  (1974)  and  Baldwin-Lomax  (1978)  zero-equation  models 
are  considered  as  standards  from  which  other  turbulence  models  can  be  judged. 
Because  the  Cebeci-Smith  model  embodies  significant  physics  for  wall  bounded 
parallel  shear  flows  that  include  pressure  gradient,  mass  transfer,  and  compressibility, 
details  of  this  model  will  be  given.  Modifications  to  include  rough  wall  effects  based  on 
the' Work  of  Laganelli  (1-75)  and  independently  discovered  (for  incompressible  flow) 
by  Krogstad  (6-91)  will  also  be  addressed.  The  reason  for  this  development  is  to 
incorporate  techniques  suggested  by  Cebeci-Smith  (1974)  and  Laganelli  (1-75  and 
12-89)  together  with  the  scales  of  velocity  and  length  as  provided  in  more  recent  work 
which  are  based  on  the  Prandtl-Kolmogorov  concept.  These  techniques  will 
encompass  wall  conditions  for  curvature,  mass  transfer,  and  roughness. 
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2.1  Zero>Equation  Turbulence  Models 


2.1.1  Cebeci-Smith  (C-S)  Model 

The  Cebeci-Smith  turbulence  model  is  a  two-layer  model  consisting  of  the 
Prandtl  mixing  length  (inner-layer)  and  the  Clauser  model  (for  the  outer-layer).  The 
inner  region  can  be  expressed  as 

-p  uV  =  p  i2  (du/dy)2  =  p  Vmi  (du/9y) 

(D 

where 

Vmi=  ilawayl  (2) 

Introducing  the  dimensionless  groups  u+  =  u/u^  y+  =  uxy/uw-  t+  =  ux  iA)w»  and 
A+=ux  A/uw*  together  with  the  van  Driest  damping  function,  the  eddy  viscosity  becomes 

^mi  =  k2(y+)2[l  -  exp  ( -y+/A+)2]  I5u+/9y+l  Vw 

The  outer  region  is  expressed  as 

Vnv,  =  0.0168  8*Uey0 

where  y0*  the  Klebanoff  intermittency  factor,  is  defined  as 

y0  =  [  1+5.5  (y/S^r1 

The  strain  gradient  I3u/9yl  can  be  expressed  as  I3u/9y  +  Sv/9xl  for  two-dimensional 

problems.  The  non-dimensional  parameter  A+  can  be  a  function  of  pressure  gradient 
as  well  as  wall  conditions  (surface  roughness  and  mass  transfer).  For  smooth  walls 

without  pressure  gradients,  the  parameter  has  been  determined  to  be  A+  =  26. 

The  van  Driest  damping  functions  is  based  on  Stokes  flow  for  an  infinite  flat  plate 
undergoing  a  simple  harmonic  oscillation  where  the  amplitude  of  motion  diminishes 
from  the  wail  by  the  factor  e*y/A.  The  constant  A  depends  upon  the  frequency  of 
oscillation  and  the  fluid  viscosity.  By  fixing  the  plate  relative  to  the  fluid  (which  now 


(4) 

(5) 
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oscillates),  a  factor  {1  -  e*y/A)  is  applied  to  the  fluid  oscillation  to  account  for  wall 
damping  effects. 

In  order  to  accomodate  compressiblity,  curviture,  and  mass  transfer  effects,  C-S 
modified  the  damping  coefficient  A+such  that 

A*  =  26  Jj*  (pw/?)1/2  (5M  (6) 

^compressible  contribution 
-^pressure  gradient  &  mass  transfer  contribution 
-^incompressible  value 

Details  of  this  development  are  provided  in  Appendix  A.  In  equation  (6),  the 

average  properties  ~p  and  p  are  usually  chosen  to  be  local  values  of  p  and  p.  the 
compressible  term  has  been  shown  to  be  a  valuable  contribution  for  treating  non- 
adiabatic  wail  conditions. 

Further  modifications  have  been  made  to  the  C-S  model  to  incorporate  wall 
roughness  (see  Appendix  A)  by  introducing  a  vortex  generating  factor  to  the  mixing 
length  as  suggested  by  van  Driest;  namely 

i£  =  Ky*  [  1  -  exp  (-yV26)  +  exp  (-6 Oy+/k?)]  (7) 

However,  the  above  does  not  allow  reoughness  heights  to  grow  for  kr+  >  60  where  the 

viscous  sublayer  would  be  eradicated  and  is  in  conflict  with  experimental  data.  In 
order  to  allow  the  damping  factor  to  be  in  excess  of  unity,  an  incrumental  factor  due  to 
roughness  is  introduced  such  that 

Dkr  =  [l+krV 30y+]  exp (-2.3^) 
and  the  damping  becomes 

D  ~  Djn»ooth+  ^kr 


(8) 

(9) 


where  D^r  is  the  sole  means  of  introducing  roughness  effects. 


The  turbulent  viscosity  coefficients  is  transitioned  from  the  inner  to  the  outer 
layer  about  a  switch  yc  point  defined  as  follows 


n,-  /  Htj@y<yc 

lpto@y>yc 


(10) 


where  yc  is  the  smallest  value  of  y  where  tyj  =  ty0. 
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2.2  Characteristic  Scales  of  Turbulence  Modeling 

In  the  previous  section,  it  was  noted  that  zero-equation  models  provided  length 
and  velocity  scales  using  prescribed  algebraic  discriptions.  One-equation  models 
used  the  turbulent  kinetic  energy  field  equation  for  the  velocity  scale  (k1/2)  and  an 
algebraic  prescription  for  the  length  scale.  Two-equation  turbulence  models  required 
field  equations  for  both  the  velocity  (k1/2)  and  length  (k3/2/e)  scales.  Here,  the  length 
scale  has  the  characteristics  required  by  Kolmogorov  (i.e.  i  =  k1/2/f  *  \&2fe  where  f  =  to 
=  e/k)  and  e  is  determined  by  the  dissipation  rate  equation 

Vt-C^tf/e  (11) 

In  the  above  equation,  the  damping  term  f^  is  critical  for  mean  flow  predictions  and 
must  have  the  correct  asympototic  dependence  as  y  -->0  in  order  for  vt  ~  0  (y3). 

since  k3/2/e  is  a  characteristic  length  scale  of  energy  containing  eddies  -0  (i) 
and  the  length  scale  k3/2/e  -  0  (y3),  a  modification  to  the  eddy  viscosity  is  required  to 
allow  for  k3/2/e~  0  (y).  For  this  situation,  the  damping  term  f^  must  be  on  the  order  of  y 

for  vt  to  be  correct  and  approach  unity  away  from  the  wall.  The  scaling  requirements 
are  as  follows  k  _  0  e  -  o  (1),  x  ~  0  (y2),  %  ~  0  (y) 

u'u'-Oty2),  u'v'~0(y3),— -0(1) 

ay  \ '*■) 

then  vt  -  0  (y3) 

Appendix  H  provides  a  further  discussion  of  scaling  issues. 

An  examination  of  the  C-S  model  shows  that 

-p  u' v' =  p  vt  — =  p  (Ky  D)2 o  (y4) 

dy  By  3y 

where  as  the  theory  of  kinematics  shows  -p  u'  v'  -0  (y3)  and  could  be  -0  (y4)  if  all 
fluctuating  quantity  correlations  were  used.  Here,  we  note  that  the  C-S  model  would 

yield  -p  u'  v'  ~0  (y2)  without  damping  and  would  provide  an  incorrect  asymptotic 
behavior  at  the  wall. 

The  Baldwin-Lomax  (B-L)  model,  see  Appendix  B  for  details,  has  similar 
features  that  are  more  avantages.  Specifically,  the  B-L  model  eliminates  the  outer 
length  scale  requirement  of  boundary  layer  thickness  and  uses  the  magnitude  of  the 
vorticity  vectoras  opposed  to  shear  strain  for  velocity  scale. 
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The  diagram  below  depicts  the  characteristics  of  the  C-S  and  B-L  models.  For 
attached  boundary  layer  flow,  both  C-S  and  B-L  models  provide  adequate  solutions  to 
compressible  flow  problems  (modest  curvature).  The  B-L  model  is  shown  to  have  a 
maximum  wake  function  (y  ©  D)  at  y/8  of  approximately  60%  where  the  wake  factor 
reaches  a  maximum  (representative  of  the  velocity  scale).  The  turbulent  viscosity  is 
shown  for  the  inner  and  outer  regions  with  the  switch  point  at  y^.  When  flow 
separation  is  experienced,  the  wake  provides  a  double  “hump”  where  the  first  hump 

exceeds  the  wake  function  maximum  value  at  y/8  of  0.6.  Here,  the  turbulent  viscosity 
can  yield  values  on  order  of  magnitude  less  at  the  switch  point  than  values  attained  for 
attached  flow  conditions. 


Both  C-S  and  B-L  models  have  been  shown  to  be  deficient  for  separated  flows. 
This  is  a  consequence  of  strong  pressure  gardients  where  separation  can  cause  the 
shear  stress  to  change  sign.  Moreover,  use  of  law  of  the  wall  relations  and  turbulence 
equations  based  on  wall  shear  stress  fail  since  the  near-waii  region  is  not  controlled 

by  u%.  Finally,  the  constants  CCp  and  Ckleb  'n  Table  B.1  have  been  shown  to  be  a 
function  of  compressibility,  wall  temperature,  and  skin-friction. 
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2.3  One-Half  Equation  Model:  Johnson  and  King  (J-K) 

The  Johnson  and  King  (11-85)  model  is  a  non-equiiibrium  eddy  viscosity  model 
that  is  formulated  with  an  algebraic  eddy-viscosity  relation  and  an  ordinary  differential 
equation  (ODE)  based  on  a  simplified  turbulent  kinetic  energy  field  equation.  This 
equation  is  based  on  the  distribution  of  maximum  turbulent  shear  stress  which  controls 
the  outer  layer  and  ensures  lag  effects  between  stress/strain  terms.  Appendix  C 
provides  details  for  this  model. 

The  inner  and  outer-layer  eddy  viscosity  relations  are  expressed  as 


Vti  =  K  y  Dj.k  (u'  v'  m)1/2 

(13) 

vto  =  a  (x)  (0.0168)  ue5  Yo 

(14) 

The  term  a  (x)  represents  an  unknown  modeling  parameter  (history  effects) 
which  is  solved  by  the  ODE,  Eq  (9)  in  Table  C.1 ,  and  is  adjusted  to  satisfy  the  relation 


(u'v'm)‘« - 

v,=  ■— —  (u'v'n,)1'2 


(du/@y) 


01 


(15) 


which  can  also  be  expressed  as 


v>  =  CtmU/2  <16> 

where  the  constant  C  and  the  length  scale  im  will  be  addressed  for  comparitive 

purposes  relative  to  one-equation  and  two-equation  turbulence  models.  Further 
discussion  concerning  the  J-K  model  and  will  be  provided. 

2.4  Characteristics  of  One-Equation  Models 

One-equation  turbulence  models  were  developed  in  order  to  incorporate  histroy 
effects  of  the  turbulent  structure.  This  is  a  consequence  of  the  inability  of  the  turbulent 
viscosity  to  follow  sudden  changes  in  the  flow  structure  such  as  shock/boundary  layer 
interactions  and  transition  from  smooth  walls  to  walls  featuring  roughness  or  blowing. 
The  eddy  viscosity  is  prescribed  as  =  vt  (k,i)  where  the  turbulent  kinetic  energy  k  is 
solved  with  a  transport  equation.  A  closure  relationship  for  the  dissipation  rate  is 
required  such  that  e  =  e  (k,  i,  v)  where  the  length  scale  i  is  prescribed 

The  inner  and  outer-region  dissipation  rate  ate  given  as 

£  ~  v  k/i2  (17) 
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and 


e  ~  k3/2/ 1 

(18) 

where  a  composite  of  two  regions  gives 

e  =  (C2k3/2/i)(l+C3RT) 

(19) 

for 

Rt  =  k1/2  i/v 

(20) 

In  the  above,  as  y  -->  0,  RT  -->  0,  and  RT  <  <  C3  giving  e  =  C2  C3  v  k/i2.  On  the  other 
hand  as  y  -->  rt  >  >  i  and  e  =  C2  k3/2/i. 

The  eddy  viscosity  is  expressed  as 

vt  =  Q  k1/2 1 D 

where  the  damping  term  D  is  expressed  as 

D  =  1  -  exp  [-C5  k1/2  y/v] 

and  the  length  scale  is  given  as 

v=  f  Ki@y^yc 
*•  C65@  y  >  yc 

A  number  of  variations  have  appeared  for  the  definition  of  of  the  length  scale  in 
the  inner  and  outer  regions  of  the  boundary  layers.  These  variation  are  presented  in 
Appendix  D  for  the  Norris  and  Reynolds  (1975)  model,  Appendix  E  for  the  Goldberg 
(6-91)  model,  and  Appendix  F  for  the  Baldwin  and  Barth  (1-91)  model. 

2.5  Two-Equation  Turbulence  Models 

The  two-equation  transport  model  for  turbulent  kinetic  energy  (k)  and  its 

dissipation  rate  (e)  have  been  formulated  to  define  appropriate  velocity  and  length 
scales  that  zero-equation  models  must  specify  in  an  ad-hoc  manner.  These  equations 
therefore  provide  a  bases  for  examining  more  complex  flow  fields.  It  is  beyond  the 
scope  of  this  paper  to  give  details  of  the  field  equation  development;  accordingly,  only 
a  discussion  of  the  results  of  the  models  will  be  provided. 

The  conservation  equations  for  a  turbulent  compressible  flow  field  are  derived  by 
expressing  the  dependent  variables  as  a  sum  of  mean  and  fluctuating  quantities.  The 
results  are  ensemble  averaged  with  the  N-S  equations  (with  variable  fluid  properties) 


(21) 

(22) 

(23) 
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together  with  a  total  energy  equation.  The  instantaneous  velocities  and  temperatures 
can  be  mass-weighted  averaged  using  the  Favre  technique  (see  Cebeci  and  Smith, 
1974)  which  provides  terms  similar  to  the  incompressible  counterpart.  For  the 
compressible  situation,  density  and  pressure  are  ensemble  averaged  without 
weighting. 

The  momentum  equation  can  be  multiplied  by  the  instantaneous  velocity  which 
is  subsequently  ensemble  averaged  to  provide  an  equation  for  the  kinetic  energy  of 

turbulent  fluctuations.  The  dissipation  rate  of  turbulent  kinetic  energy  (e)  can  be 

expressed  in  analogy  to  the  turbulent  kinetic  energy  format  that  provides  the  k-e 
turbulent  transport  equations.  These  field  equations  which  are  provided  in  a  number 
of  references  stated  above  will  be  expressed  in  terms  of  mean  flow  quantities.  The 
Appendix  (G  and  H)  of  this  report  provides  the  conservation  equations,  transport 

equations  and  constitutive  relations  for  the  k-e  two-equation  model  (Wilcox,  6-91)  as 
well  as  the  k-e  two-equation  model  (Speziale  et  al,  6-90).  The  latter  retains  the 
fluctuating  quantities  for  the  purpose  of  discussion  for  embodied  physics.  The 
resulting  field  equations  are: 

V.(v+^)  V]c  +  p-e  =  Ek  (24) 


where 


DE  =  v.(v4Vl)  v£  +Qf Jp  -c2f&-  +  Ee  +  E 
Dt  Og  K  K 

D  (  )/Dt  =  3  (  )/at  +  U*  v 

vt  =  C(if^k2/e 


/9u,  9u;\3uj  2  /^ukY 


e  =  e  +  DE;  e  =  v 


d  u'i  d  u'j 
— — ■ 

d  Xj  d  Xj 


Ek  =  CD  f  (Mt)  P  e 


(25) 


(26) 

(27) 

(28) 


EE  =  f(Mt) 


The  functional  form  of  the  turbulent  Mach  number  f(Mt)  has  several  prescriptions, 
(Dash,  6-91),  (Sarkar  et  al,  11-89),  (Zeman,2-90),  and  (Wilcox,  6-91). 
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The  standard  high  Reynolds  number  form  of  the  k-e  field  equations  provides  for 
the  following  definitions: 

The  constants 


=  1 .0,  a£  =  1.3,  Ci  =  1.44,  C2  =  1.92,  =  0.09 

The  functions 

f1  =  f 2  =  1^*1,  represent  wall  functions  for  low  Re  flows 
Ek  =  Ee  =  O,  represents  compressibility  corrections 

E  =  O,  source  term  to  dissipation 
De  =  O,  correction  to  dissipation 


(29) 


The  eddy  viscosity  is  given  by  equation  (11). 


As  previously  noted  the  two  most  popular  two-equation  turbulence  models 
consist  of  the  Jones  and  Launder  (2-72)  k-  e  model  and  the  Saffman  (1970)  and 

Wilcox  (1-91  and  6-91)  k-to  type.  The  common  field  equation  of  the  two-equation 
models  is  the  turbulent  kinetic  energy  (TKE)  equation,  eq.  (24).  The  rate  of  dissipation 

of  TKE  (e)  is  determined  from  field  equations  that  are  similar  to  the  TKE  transport 
equation.  The  relationship  of  the  velocity  and  length  scales  with  eddy  viscosity  are 
given  as 


Pt  =  Cn  fyp  k1/2 1  =  Cji  p  k2/e  =  fn  p  k/co 


(30) 


for 

i  =  k1/2/  a) ,  a  =  e/k 

where  f^  and  have  been  previously  defined  and  to  is  the  specific  dissipation  rate. 

The  standard  k-  e  model  has  been  used  to  predict  a  wide  variety  of  flows  that 
include  jets,  wakes,  mixing  layers,  boundary  layer  (attached  and  separated),  and  flows 
experiencing  wall  mass  transfer  and  roughness.  Due  to  its  inherent  limitation  with  the 

e  transport  equation  behavior  approaching  solid  boundaries,  wail  techniques  have 
been  required  to  provide  solutions  to  a  wide  range  of  engineering  problems.  Since 
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ad-hoc  techniques  are  required  for  the  near-wall  region  which  can  be  numerically  stiff 
in  turbulent  boundary  layer  flows,  alternate  forms  of  two-equation  models  have  been 
sought.  These  include  the  k-oo  model  of  Wilcox  (1-91),  the  k-W  model  of  Spalding  (11- 
82)  and  the  k-t  model  of  Speziale  et  al  (6-90).  The  latter  two  models  are  essentially 

variations  of  the  Wilcox  k-co  model.  Appendix  G  and  H  of  this  report  provide  the  field 
equations,  constitutive  relations,  conservation  of  mass,  momentum  and  energy  as  well 
as  closure  coefficients  for  these  models. 

As  noted  in  the  introduction,  a  position  on  choice  of  turbulence  closure 
modeling  is  required  in  order  to  meet  near-term  goals  of  providing  design  resolution  to 
practical  engineering  problems.  Turbulence  modeling  at  the  two-equation  level  meets 
this  requirement.  The  k-  e  two-equation  model  has  been  clearly  the  most  used  and 
has  developed  an  extensive  experience  base  to  meet  the  challenges  of  many  complex 
flow  problems.  This  experience  base  also  includes  the  adaptation  of  the  k-  e  models 
into  multizone  PNS  solvers  which  has  not  been  demonstrated  by  the  k-o)  type. 
Moreover,  coupling  of  unit  level  problems  has  been  demonstrated  by  the  k-  e  type  with 

proper  use  of  wall  techniques,  whereas  the  k-co  type  require  boundary  condition 
adjustments  (usually  with  log-law  relations  to  account  for  wail  effects  (roughness  and 
blowing  for  example)  which  can  present  problems  under  compressible  flow  conditions. 

With  the  recent  success  of  the  two-layer  models  that  incorporate  the  standard  k- 
e  model  and  one-equation  model  for  length  scale  in  the  near-wall  region  thereby 
eliminating  difficulty  of  the  e  equation  in  approaching  solid  wall  boundaries,  it  is  our 
position  that  the  k-e  two-equation  type  model  can  meet  near-term  goals. 

It  will  be  demonstrated  that  a  hybrid  model  that  features  the  standard  k-e 
field  equations  for  the  outer  layer  and  a  one-equation  model  for  the  inner  near-wall 
region  provides  the  correct  asymptotic  behavior  for  dissipation  and  eliminates 
numerical  stiffness  problems.  The  experience  base  with  use  of  other  wail  techniques 
as  well  as  adapting  corrections  for  compressibility  and  curvature  provide  the  capability 
to  treat  a  wide  class  of  engineering  problems  that  include  free  shear  layer,  mixing 
layers,  and  wall  bounded  flows  including  separation  subject  to  compressible  flow 
conditions. 


Table  II  provides  a  comparison  of  the  J-K  model  together  with  the  C-S  model  and 
the  N-R  model.  Since  the  Rodi  hybrid  model  is  based  on  the  N-R  one-equation  model 
for  the  length  scale  in  the  near  wall  region,  it  has  been  added  for  comparison.  The 
inner-layer  eddy  viscosity  has  been  structured  in  terms  of  fluctuating  components-  in 
order  to  make  a  one  to  one  comparison  of  ail  models.  Three  distinctive  variations 
occur  in  length  and  velocity  scales  that  are  the  underlying  reasons  why  the  J-K  and 
Rodi  hybrid  models  have  had  demonstrated  success  in  both  incompressible  and 
compressible  flow.  The  most  important  is  the  choice  of  velocity  scale.  The  C-S  model 
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Table  II.  COMPARISON  OF  TURBULENCE  MODEL  LENGTH/VELOCITY  SCALES 
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ons  in  the  switch  point  locations  are  prescribed  in  the  models 


uses  the  wall  shear  stress  velocity  (u  )  while  the  J-K  and  N-  R  use  a  variation  on 
the  Prandtl  concept.  The  J-K  model  uses  (-u'vj172  while  the  N-R  uses  (k)1/2  where  k  is 
the  sum  of  the  square  of  the  stress-strain  diagonal  components  i.e.1/2  [(u1)2  +  (v')2  + 

(w)2].  A  variation  in  the  damping  coefficient  A+  also  exists  as  well  as  the  exponent  of 
the  damping  term  (i.e.,  J-K  uses  a  power  of  2,  while  the  other  features  a  power  of 
unity).  This  was  required  to  allow  for  damping  to  adjust  to  -0  (y)  instead  of  -  0  (y2), 
note  the  power  of  unity  on  the  length  scale  icy.  With  reference  to  eqs  (15)  and  (16),  a 
comparison  between  the  length  scales  of  the  N-R  and  J-K  models  can  be  made.  Here, 
im  compares  directly  with  the  N-R  i^,  eq  D-2,  with  the  velocity  scales  as  described 
above.  Variations  in  the  constants  are  also  noted.  It  would  be  interesting  to  re¬ 
examine  the  J-K  model  with  the  correct  value  for  the  damping  coefficient  A+  (=26)  and 
allow  for  (k  y  D)2  as  developed  in  the  zero-equation  and  one  equation  models. 

The  elegance  of  the  J-K  model  is  provided  in  the  combination  of  the  use  of  a 

turbulence  energy  term  for  the  velocity  scale  in  the  inner-layer  and  the  parameter  a  (x) 
which  is  determined  from  an  ODE  for  maximum  Reynolds  shear  stress  in  the  outer- 
layer.  The  outer  layer  solution  uses  an  ODE  while  the  N-R  model  uses  a  turbulent 
kinetic  energy  field  equation  with  a  modified  Kolmogorov  term  for  dissipation.  Since 
the  characteristic  scales  are  based  on  maximum  shear  conditions,  wall  effects  such  as 
pressure  heat-  transfer,  and  drag  are  reasonably  predicted.  However,  details  of  the 
shear  layer  structure  cannot  be  obtained  by  the  J-K  model  as  would  be  provided  by 
the  N-R  model  and  to  a  greater  extent  by  the  Rodi  hybrid  model. 

2.6  Higher  Order  Steady  State  Turbulence  Models 

Beyond  the  two-equation  level,  a  broad  spectrum  of  turbulence  model  have 
been  formulated  to  remedy  fundamental  deficiencies  such  as  Boussinesq  gradient 

transport  hypothesis  -pu‘  x>'  =  —  and  the  assumption  of  isotropy.  The  simplest 

x  3y 

remedy  for  removing  the  isotropy  assumption  can  be  performed  at  a  two-equation 
level  with  significant  additional  algebra  and  is  called  Algebraic  Reynolds  Stress  (ARS) 
model.  Such  models  were  first  formulated  in  the  mid-to-later  seventies  at  Imperial 
College  and  have  been  predominately  employed  for  vortical  turbulent  flows  such  as 
corner  flow  and  V/STOL  jets.  Until  quite  recently,  all  ARS  work  has  emphasized 
incompressible  flows.  In  the  recent  work  of  Burr  and  Dutton  (1-90)  at  the  University  of 
Illinois,  ARS  models  have  been  extended  to  highly  compressible  flows  with  a  primary 
emphasis  on  NASP-oriented  flow  problems. 

The  next  class  of  higher-order  turbulence  models  are  multi-scale  models  which 
distinguish  between  turbulent  producing  large  eddies  and  turbulent  dissipating  small 
eddies.  These  are  four  equation  models  which  solve  for  the  turbulent  kinetic  energy  of 
the  large  and  small  eddies  and  turbulent  dissipation  rates  for  the  transfer  from  large  to 
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small  scale  in  addition  to  the  basic  dissipation  rate.  Several  models  have  been 

formulated  as  extensions  for  the  k-e  model  with  pioneering  work  in  this  area  again 
performed  at  Imperial  College  (Hanjalic  and  Launder)  and  more  recently  by  Kim  and 
Chen  (1-88)  at  the  University  of  Alabama/Huntsville.  This  work  too  has  focused  on 
incompressible  flows.  Attempts  to  deal  with  compressible  flows  are  provided  by  Adol- 
Hamid  and  Wilmoth  (6-87);NASA/LRC,  using  the  compressibility-correction  to  the  k-e 
model  of  Dash  et.  al  (8-75)'  Section  3.3  addresses  compressible  issues. 

Multi  scale  work  with  the  k-oo  model  has  also  been  performed  as  described  by 
Wilcox  in  several  recent  AIAA  Journal  articles  (11-88  and  7-90).  Multiscale  models 
serve  to  eliminate  the  Boussinesq  approximation  since  they  compute  the  Reynolds 
Stress  directly. 

Next  in  the  hierarchy  of  higher  order  models  are  the  Reynolds-Stress  models 
which  solve  partial  differential  equations  for  u‘  u‘ ,  v*  v* ,  w*  w\  and  u'  v'  in  conjunction 
with  a  dissipation  rate  equation.  These  models  have  focussed  primarily  on 
incompressible  flows  until  quite  recently.  The  new  models  of  Sarker  et.  al  (6-90)  at 
NASA/LARC  and  Zeman  (2-90)  at  Stanford  are  compressible  extensions  of 
conventional  Reynolds-Stress  models  which  use  Favre  (density-weighted)  averaging 
to  account  for  density  fluctuations.  The  basic  compressible  models  do  not  account  for 
observed  compressibility  effects  and  must  be  supplemented  by  modifications  to  the 
dissipation  rate  equation  (compressible-dissipation  extension)  to  properly  account  for 
observed  compressibility  effects. 

The  most  sophisticated  of  the  steady  state  turbulence  models  are  full  second- 
order  closure  models  which  solve  partial  equations  (pde’s)  for  all  second-order 
correlated  terms  such  as  pur,  .  .  .  ,  h'  u\  ....  etc.  Third-order  correlations  are 
expressed  as  sums  of  the  2nd-order  correlations.  The  leading  proponent  of  this  work 
has  been  Donaldson  (1966)  at  ARAR  Early  issues  involved  the  use  of  algebraic 
length  scale  equation  which  has  since  been  remedied.  Such  models  lost  momentum 
in  conclusions  drawn  in  1982  AFOSR/Stanford  Turbulence  Model  Conference  where 
the  claimed  “invariance”  of  coefficients  was  questioned  and  refuted.  For  flows  with 
multi-component  species,  the  number  of  pde’s  can  be  overwhelming  making  such 
methodology  highly  impractical. 
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3.0 


TURBULENCE  MODEL  RECOMMENDATIONS 


3.1  Standard  High  Reynolds  Number  k-e  Model 

As  previously  noted,  our  position  on  choice  of  a  turbulence  closure  model  is  the 
standard  two-equation  k-e  model  and  its  subsequent  modifications  to  treat  a  wide 
range  of  engineering  problems.  The  extensive  experience  base  generated  by  these 
models  should  provide  a  level  of  generality  with  proper  definitions  for  closure 
constants  (or  parameters)  together  with  near-wall  terms  to  handle  multiple  flow  states. 
Moreover,  the  versatility  of  the  models  and  adaptation  to  existing  industrial  codes 
allows  for  robust  transfer  of  the  methodology. 

Two-equation  models  also  provide  other  obvious  features  over  lower-level 
models.  This  includes  the  use  of  better  physics  in  the  field  equations  to  handle  the 
outer  layer  region  where  accelerated  flows,  severe  pressure  gradients,  and  mass 
transfer  (diffusion)  can  have  an  impact.  Details  of  the  turbulence  structure  are 
provided  which  can  help  provide  insight  for  a  wall  function  modeling.  The  models  are 
not  computationally  intensive  and  are  readily  adaptable  to  current  industrial  codes. 

Equations  (24)  through  (29)  presented  the  generalized  two-equation  model  with 
and  without  modifications.  Table  III  has  been  structured  to  show  how  these 
modifications  have  been  made  to  the  standard  high  Reynolds  number  equations. 
These  modifications  include  near  wall  techniques,  compressible  corrections,  and 
geometric  (curvature)  corrections.  Each  of  these  modifications  will  be  discussed.  The 

closure  constants  that  are  recommended  for  the  standard  k-e  model  were  acquired 
from  incompressible  experiments  of  free  turbulent  flows  that  include  free  shear  layers 
(jets)  as  well  as  confined  mixing  layers.  These  same  constants  have  been 
recommended  for  boundary  layer  flows.  The  relative  merit  of  these  constants  in  flows 
subject  to  wall  bounded  flows  which  has  necessitated  wall  damping  functions  requires 
further  investigation  relative  to  compressible  effects.. 

3.2  Near-Wall  Techniques 

Four  levels  of  near-wall  techniques  have  been  characterized  as  wall  functions, 
bridging  techniques  using  van  Driest  damping,  low  Reynolds  number  techniques,  and 
hybrid  techniques.  The  hybrid  techniques  actually  feature  two-layer  models  that 

incorporate  the  standard  k-e  model  in  the  outer  layer  with  a  one-equation  model  for  the 
inner  layer.  Wall  functions  use  law  of  the  wall  equations  to  match  a  grid  point  usually 
taken  at  the  edge  of  the  viscous  sublayer  with  no  attempt  to  resolve  the  sublayer 

region.  This  point  is  then  used  as  a  boundary  condition  for  the  outer  region  for  the  k-e 
equations.  Modified  versions  of  the  k-equation  are  also  used  at  this  point  to  solve  for 
the  kinetic  energy  as  opposed  to  setting  its  value.  Wall  shear  stress  velocity  is  taken 
as  the  characteristic  velocity.  These  techniques  have  the  advantage  of  not 
encountering  numerical  problems  (no  grid  points)  in  the  sublayer  and  can  save 
computational  time.  However,  they  require  ad-hoc  techniques  between  the  wall  and 
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Table  HI.  TWO-EQUATION  k-e  MODEL  VARIATIONS 


-  'Mv-7p  *  Ci',P  -  *  E  ♦  E, 

V,  -  C^k »/€ 

e  ®  c  +  DE 


STANDARD  HIGH  Re  k-e  MODEL 

0|f  —  If  o(  —  1«3|C1  —  1.44,  C, 

=  1.92,  &  C,  =  0.09 

fi  =  f2  -  f»,  =  unity, 

Low  Re  Wall  Functions 

E*  =  E,  -  0, 

Compressible  Corrections 

o 

m 

II 

O 

** 

Correction  to  Dissipation 

E  =0, 

Source  Term 

DASH 

SARKAR  et  at 


WALL 

FUNCTIONS. 

-  Spaldlng-Laundor 

-  Chian  g-Laundar 

-  Viogas-Rubasin 

-  Gondd 


NEAR  WALL 
TECHNIQUES 


GEOMETRIC 

CORRECTIONS 


-  AXISYMMETRIC 

-  STREAMWISE  CURVUURE 

-  SWIRL 

-  VORTICAL 


BRIDGING 

TECHNIQUES 

(Yan-Diteat) 
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•  Aroraotal 
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•ECHNIQUES 

•  Jonss-Laundar 

-  Chian 
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TECHNIQUES 
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(Van  Driost) 
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grid  first  point  which  do  not  have  established  physics.  Moreover,  use  of  law  of  the  wail 
relations  has  been  demonstrated  to  be  inadequate  for  complex  flows  where 
inaccuracies  affect  turbulence  quantities  as  well  as  mean  flow  properties. 

A  second  approach  uses  bridging  techniques  between  the  wall  and  grid  point  at 
the  edge  of  the  viscous  sublayer  with  consideration  to  multiple  layers  in  between. 
Here,  the  van  Driest  damping  functions  are  applied  to  provide  algebraic  relations  to 
solve  for  the  flow  conditions.  While  this  method  provides  better  physics,  the 
phenomena  is  linear  and  does  not  account  for  lag  between  the  stress  and  mean  flow 
properties.  An  advantage  of  this  technique  however,  is  the  ability  to  incorporate  wall 
conditions  such  as  roughness,  blowing,  and  pressure  gradients  through  modifications 
of  the  van  Driest  damping  function  (or  coefficient).  This  has  been  demonstrated  by 
Aurora  et  al  (  11-82)  as  well  as  Cebeci-Smith  (see  Appendix  A). 

The  third  technique  considers  use  of  the  field  equations  extended  to  the  wall 
with  low  Reynolds  number  techniques  that  require  damping  function  adjustments  to 
the  constants,  as  well  as  to  the  dissipation  term.  In  some  instances,  the  addition  of  a 
source  term  to  the  dissipation  has  been  used.  These  techniques  provide  better 
physics  and  have  been  shown  to  predict  transitional  boundary  layer  flow  as  well  as  re- 
laminarization.  In  some  cases,  a  non-physical  form  of  dissipation  is  solved  within  the 
viscous  sublayer  requiring  a  modified  kinetic  energy  equation  to  balance  the 
dissipation.  As  in  the  van  Driest  technique,  exponential  damping  is  used  to 
compensate  errors  for  the  new  ad-hoc  technique  which  do  not  allow  for  the  eddy 

viscosity  to  correctly  asymptote  to  the  wall  (i.e.  it  does  not  have  an  O  (y2)  behavior). 
These  techniques  also  offer  opportunities  to  investigate  wail  conditions  of  roughness, 
blowing, and  curvature  through  the  damping  terms  or  damping  coefficients.  One 
interesting  aspect  of  the  low  Reynolds  number  methods  is  the  choice  of  velocity  scales 
which  use  both  the  wall  shear  velocity  as  well  as  a  velocity  based  on  kinetic  energy. 

The  fourth  technique  considers  a  hybrid  method  that  uses  the  standard  k-e 
model  in  the  outer  layer  of  a  two-layer  model  and  a  one-equation  model  for  the  inner 
layer.  This  approach  has  shown  recent  success  for  both  incompressible  and 
compressible  flows  to  treat  severe  pressure  gradient  regions.  Moreover,  the  one- 
equation  model  eliminates  the  difficulty  of  the  limitation  of  the  dissipation  to  approach 
solid  boundaries  and  provides  for  techniques  to  treat  roughness,  mass  transfer, 
curvature,  and  transition.  The  technique  used  by  Rodi  (1-91)  for  incompressible  flow 
and  extended  to  compressible  flow  by  Horstman  (6-91)  uses  the  one-equation  model 
of  Norris-Reynolds  (see  Appendix  D)  and  will  be  the  candidate  methodology  for  a 
generalized  turbulence  model. 


3.2.1  RODI  Hybrid  Model 

The  Rodi  model  (1-91)  combines  the  k-e  transport  equations  with  a  one-equation 
length  scale  model  for  the  near  wall  region.  The  rationale  and  features  of  this  model 
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have  been  discussed.  Table  IV  provides  the  model  characteristics.  The  bulk  of  the 

shear  layer  is  determined  by  the  k-e  equations.  As  evidenced  by  this  review,  a  number 
of  variations  have  been  used  for  the  damping  functions  for  the  one-equation,  length 
scale  as  well  as  the  matching  criteria  for  the  inner/outer  layers,  it  has  been 
demonstrated  that  the  use  of  the  N-R  one-equation  model  provides  results  for 
incompressible  flow  that  are  consistent  with  low  Reynolds  k-e  models,  and  have  been 
shown  to  be  better  for  flows  with  adverse  pressure  gradients.  Moreover,  the 
recent  work  of  Horstman  (6-91)  has  shown  that  the  Rodi  hybrid  model  provides 
resolution  to  compressible  2D/3D  shock/boundary  layer  interactions  for  supersonic 
and  hypersonic  conditions. 

The  standard  high  Reynolds  number  k-e  model  is  used  for  the  bulk  of  the  shear 
region  away  from  the  wail  and  calculates  the  distribution  of  eddy  viscosity  from  (see 
eq.  (11)). 

Dt  =  C^k2/e 

The  near  wall  region  is  formulated  by  a  one-equation  model  whose  length  scale 
is  prescribed  by  Equations  (D-1)  through  (D-4).  It  should  be  noted  that  the  length 
scales  and  i£  are  the  same  in  the  log-law  region  when  the  coefficient  is  chosen 

as  the  square  of  u'v'/k  under  equilibrium  conditions.  In  the  near  wall  regions, 
gradients  require  that  i„  and  i£  be  different  and  a  reduction  in  y  can  be  tailored  after 

the  van  Driest  formulation  with  a  proper  definition  of  velocity  scales.  The  N-R  model 
provides  for  such  a  function. 

Several  other  phenomena  make  the  Rodi  model  attractive  for  potential 
universality.  This  includes  the  ability  to  provide  wall  conditions  such  as  blowing, 
roughness,  and  curvature  using  the  damping  coefficients  (C3  in  Equations  (D-2)  and 

(D-8))  as  developed  by  C-S,  Rodi  (1984),  and  Laganelli  (1-75),  as  well  as  handle 
adverse  pressure  gradient  and  compressible  flow  conditions.  In  essence,  the  two- 
layer  model  appears  capable  of  providing  solutions  of  coupled  problems.  Another 
attractive  feature  is  that  the  model  considers  a  modified  version  of  the  Kolmogorov 

function  and  provides  for  correct  length  scale  where  1  ~  y  at  the  wall.  Finally,  since  the 
transport  equations  are  being  used,  details  of  the  flow  structure  can  be  obtained  that 
are  not  available  with  lower  order  models. 


3.3  Compressible  Modifications 

The  k-e  model  was  originally,  formulated  from  a  purely  incompressible 
viewpoint.  Its  inability  to  deal  with  high  Mach  number  compressibility  effects  was  first 
made  evident  at  the  1972  NASA  Free  Shear  Flow  Conference  (NASA-72)  where  its 
performance  versus  the  classic  isoenergetic/single-stream  shear  layer  data  showed 
no  decrease  in  mixing  with  increasing  Mach  number,  whereas  the  data  indicated  a 
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Table  IV.  RODI  TWO-LAYER  (RTL)  TURBULENCE  MODEL 

(from  Horstman,  6-91) 


The  two-equation  k  -  e  turbulence  model  is 
expressed  as 

2 

- pkdivu-p€+ diffusion +  LRT  (1) 


PjSf^pVkL 

(5) 

fjt  a  I  -  exp(-0.0198  Ry) 

(6) 

Ry  =  pVk  y  /  p 

(7) 

Dpe 

Dt 


€ 

k 


-C[  — pk  div  u 


-C2f2p — + diffusion  +  LRT 


(2) 


where  y  is  the  normal  distance  to  the  wall.  The  dissipation 
rate  e  is  given  by 


e 


13.2 

pVkL/p, 


(8) 


Sy-uij  +  ujj  (3) 


Pt  ■  C^pk2/e  (4) 


where  LRT  refers  to  the  low  Reynolds  number  terms  used 
near  solid  surfaces.  The  values  used  for  the  constants 
wereCi  =  1-44.C2*  1.92,  and  Cy.= 0.09. 

The  second  model  used  here  is  the  two-layer  model 
developed  by  Rodi  and  his  co-workers  For  this  model 
the  e  equation  is  replaced  by  an  analytical  expression 
near  the  wall.  In  the  wall  region  this  model  determines  the 
eddy  viscosity  from  the  relation 


The  length  scale  L  is  assumed  proportional  to  the  dis¬ 
tance  from  the  wall,  i.e..  L  =  Cry.  The  wall  region  model 
is  matched  with  the  full  k  -  €  model  at  a  distance  from 
the  wall  where  Ry  =  150  (fjx  **  0.95).  At  this  location  the 
eddy  viscosity  given  by  Eqn.  4  is  equated  to  Eqn.  5  to 
determine  the  proportionality  constant  Cr.  For  typical 
shock-wave/boundary-layer  interaction  flows  solved  here 
Cr  varied  from  0.8  to  9.0  at  various  locations  in  the  flow 
field.  The  lowest  values  occurred  upstream  of  the  interac¬ 
tion  and  the  highest  near  rcattachmcnt  For  the  present 
computations,  maximum  and  minimum  limits  of  2.4  and 
\2  were  used  for  Cj.  With  these  limits  imposed,  the  val¬ 
ues  of  eddy  viscosity  are  no  longer  matched  between  the 
inner  and  outer  models  for  large  regions  of  the  flow  field. 
A  few  solutions  were  obtained  using  Cr  =  2.4  which  is 
equivalent  to  using  the  Prundtl  length  scale.  (Rodi 
used  a  constant  value  of  2.45  in  his  incompressible  com¬ 
putations).  Solutions  were  also  obtained  for  several  inter¬ 
action  flows  varying  the  match  point  from  Ry  =  100  to 
200.  These  results  differed  from  each  other  by  less  than 
3  percent 
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substantial  decrease  at  Mach  numbers  above  unity  (e.g.,  a  factor  of  4  decrease  at 
Mach  5).  Attempts  to  remedy  this  deficiency  has  focussed  on  free  shear  flows 

(jets/plumes/wakes)  since  the  performance  of  the  k-e  model  for  basic  boundary  layer 
flows  was  found  to  be  adequate  to  Mach  5  and  beyond  (depending  on  wall 
conditions). 

Over  the  years,  a  variety  of  compressibility-corrections  have  been  formulated  to 
address  the  issue  of  reduced  mixing  at  higher  Mach  numbers.  Explanations  of  the 
phenomenology  causing  this  behavior  have  been  varied  and  based  on  stability 
considerations,  eddy-shocklet  concepts,  etc.  The  most  widely-used  correction  for 

jets/plumes  has  been  that  of  Dash  et  al.  (8-75)  to  the  k-e  model.  The  k-e  cc  model  (i.e. 
k-e  with  compressible  corrections)  uses  a  correction  to  the  turbulent  viscosity,  eq.  (11), 
of  the  form: 

M-ti  =  Cju  K  (MT)p 

where  C^j  is  the  incompressible  coefficient  (=.09)  and  K  (Mt)  =  C^/C^j  is  a  curve  fit 
which  enforces  agreement  with  the  isoenergetic  single  stream  shear  layer  data.  In 
equation  (31),  M  is  the  turbulent  fluctuation  Mach  number 

w 

Mr  =q/a  (32) 

where  q  =  (UjUj)1/2  and  the  peak  value  of  Mx  is  used  at  each  axial  position  of  the  shear 

layer.  A  number  of  compressibility  corrections  to  the  k-e  and  higher  order  models  have 
been  proposed,  some  of  which  are  summarized  in  Table  V.  The  recent  compressible 
dissipation  corrections  of  Sarkar  et  al.  (6-90)  and  Zeman  (2-90)  appear  to  have 
generated  the  most  enthusiasm  as  evidenced  by  recent  papers  at  the  AIAA  Fluids  and 
Plasmadynamics  Meeting  in  Hawaii  (June  1991). 

Papers  by  Dash  (6-90),  Viegas  and  Rubesin  (6-91) ,  and  Wilcox  (6-91)  all  have 

explored  the  inclusion  of  these  compressible-dissipation  models  into  k-e  and  k-co 
turbulence  model  frameworks.  In  both  the  Sarkar  and  Zeman  models,  the  turbulent 

dissipation  rate,  e,  is  taken  to  be  comprised  of  a  solenoidal,  incompressible 
component,  e,.,  and  dilatational,  compressible  component,  e,,,  and  thus: 

pe  =  p(es  +  ec)  (33) 

which  replaces  ]d  e  in  the  standard  k-e  equations  (See  Table  III). 

Both  Sarkar  and  Zeman  integrate  the  incompressible  form  of  the  dissipation  equation 
to  obtain  eg  and  model  £c  as  follows: 
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TURBULENCE  MODELS  WITH  COMPRESSIBILITY-CORRECTION 
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isurr  ana  uuiton  aks  Extension  or  ke 

(1990,  AIAA  Paper  90-1463) _ 

Wilcox  kW  with  Multi-Scale  Extensions  and 

(1991,  AIAA  Paper  91-1785)  Sarkar/Zeman  Blend 


Sarkar  et  al: 


Zeman; 


where 


for  Mx>0.1,  and, 


£c  =  at  e,  M| 


£c  =  Q  F  (MT)  e$ 


F  (Mt)  =  1  -  exp  [  (-Mr  -  0.1]/  0.6)2] 


F  (Mt)  =  0 


for  M_  <  0.1.  In  their  calibrations  of  these  compressible-dissipation  models,  Sarkar 
found  (*!  =  1 .0  to  provide  the  best  agreement  and  Zeman  found  Cd  =  .75  to  work  best. 

Sarkar  and  Zeman  have  incorporated  their  compressible-dissipation  models 
into  a  Reynolds-Stress  model  framework  where  they  solve  the  field  equations  for  the 
varied  stress  components  (i.e.  for  u'u\  V7,  uV,  — ).  Simplifying  their  models  to  work 

in  a  k-e  two-equation  mode!  framework  leaves  open  some  ambiguities  with  regard  to 
the  value  of  the  turbulent  viscosity  to  be  utilized.  The  initial  interpretations  of  Dash 

(6-91)  and  Viegas/Rubesin(6-91)  differed  substantially  and  thus  produced  different 
results.  This  descrepency  was  addressed  by  Dash  et  al  (9-91 )  and  was  associated 

with  the  interpretation  of  the  turbulent  viscosity,  px,  appearing  in  the  diffusion  terms  of 

the  mean  flow  equations,  in  the  diffusion  terms  of  the  k  and  eg  equations,  and  in  the 
turbulent  production  term,  a.  Compressible  and  incompressible  values  of  p_  (j x,n  and 

X  XC 

respectively)  are  defined  as  follows: 

fi-tc  =  QioP  k2/  (Sj  +Sc)  (381 


Pts  =  C^op  k2/6s  (39) 

Dash  used  the  compressible  for  all  diffusion  terms  (mean  flow,  k  and  e),  the 
compressible  for  in  the  k  equation  and  the  incompressible  for  p  in  the  e 

X  X 

equation.  Viegas/Rubesin  used  the  incompressible  p  for  all  terms.  After  reviewing 


these  two  initial  attempts,  Dash  explored  this  matter  in  greater  detail  to  ascertain  which 
interpretation  was  most  consistent  with  the  Reynolds-Stress  formulation  of  Sarkar  et  al. 
A  matrix  of  plausible  interpretations  is  summarized  in  Table  VI. 

In  Table  VI  Df,  Dj<,  and  D£  refer  to  the  diffusive  terms  in  the  mean  flow,  k  and  e 
equations  respectively,  p„  and  p  refer  to  the  turbulent  production  term  in  the  k  and  e 

It  fc» 

equations. 

In  reviewing  the  Sarker  Reynolds-Stress  equations,  method  A  (Dash  6-91) 
appeared  inconsistent  since  different  values  of  production  were  employed  in  the  k  and 
e  equations.  Method  B  (Viegas/Rubesin,  6-91)  also  appeared  inconsistent  since 
wherever  a  gradient  transport  hypothesis  was  employed  by  Sarkar  et  al.,  a 
compressible  value  of  p  was  utilized  whereas  Viegas/Rubesin  used  incompressible 
values  throughout.  The  different  methods  for  including  compressible-dissipation  into 
k-  e  listed  in  Table  VI  were  assessed  by  Dash  via  comparisons  with  the  spread  rate 
data  correlation  of  Birch  and  Eggers  (NASA,  1 972)  LaRC  correlation,  for  asymptotic 
isoenergetic/single-stream  shear  layers.  The  Sarkar  formulation  was  used  with  cti  = 
1 .0.  Figure  1  exhibits  the  predictions  for  variants  A  -  E  vs  the  Mach  number  LaRC 
correlation.  The  spread  rate  parameter,  or  (=1.855  Ax/Ay)  is  plotted  vs  the  Mach 
number  of  the  moving  stream.  As  the  Mach  number  increases  above  sonic,  the  spread 
parameter  increased  from  the  incompressible  value  (a0  -  11)  indicating  a  marked 
decrease  in  spread  rate  with  increasing  Mach  number.  Variants  A  and  E  underpredict 
the  rate  of  mixing  throughout  (cr’s  are  too  large)  whereas  variants  B,  C  and  D 
underpredict  the  rate  of  spread  of  Mi  <1.5  and  overpredict  the  rate  of  spread  for  Mi  > 
1 .5.  Variants  C  and  D  have  identical  behavior  whereas  variant  B  shows  somewhat 
smaller  spread.  Comparing  with  Sarkar*s  Reynolds-Stress  (RS)  predictions,  variants 
A  and  E  can  be  ruled  out  whereas  variants  B,  C  and  D  ail  produce  spread  behavior 
which  appears  quite  comparable.  Version  D  has  been  selected  for  reasons  discusses 
by  Dash  (9-91). 

-  -  As  shown  in  Figure  1 ,  Variant  D  overpredicts  spread  below  M  =  1 .5  and 
underpredicts  above  M  =  1.5  (this  is  true  for  both  the  Sarkar  and  Zeman  formulations). 
In  view  of  its  simplicity,  Dash  (9-91)  chose  to  work  with  the  Sarkar  formulation  in  the 
confines  of  Variant  D,  and,  to  revise  it  to  better  match  the  spread  data.  The  revised 

version  is  described  below.  The  compressible-dissipation,  ec,  is  redefined  as: 

ecseJaMj  +  PMj]  (40) 

where  a  Zeman  lag  is  inserted  (Bf  =  M^  -X)  and  a  fourth-order  term  is  added. 

The  coefficients  utilized  are: 
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Mach  Number 


Figure  1.  Comparison  of  Variants  A-E  of  Table  VI  in 
Producing  LaRC  Isoenergetic/Single-Stream 
Spread  Sheet 


Figure  2.  Comparison  of  ke  CD  (variant  D  of  Sarkar) 
and  kcCC  Models  versus  LaRC  Spread  Sheet 
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1  (same  as  Sarkar  et  al.) 

0.1  (same  as  Zeman) 

60  (fits  LaRC  data  the  best) 

Figure  2  shows  performance  of  the  revised  model  with  the  LaRC  spread  data 
and  with  the  keCC  model  (Dash  et  ai.,  8-75),  also  calibrated  to  fit  this  same  area. 

Two  turbulence  models  have  been  calibrated  to  reproduce  the 
isoenergetic/single-stream  shear  layer  data.  We  can  designate  the  revised  Sarkar 

model  in  k-e  using  Variant  D  is  designated  as  the  keCD  model.  A  preliminary  attempt 
was  made  by  Dash  (9-91)  to  distinguish  between  their  performance  first  using  other 
data  sets,  and  then  profile  comparisons  including  the  behavior  of  turbulent  quantities. 
In  analyzing  the  isoenergetic/two-stream  data  of  Chinzei  et  al.  (1986)  and  the  non- 
isoenergetic/two-stream  data  of  Goebel  and  Dutton  (4-91),  the  spread  behavior  is 
nearly  identical  and  no  distinguishing  characteristics  are  observed.  In  comparing 
profiles,  however,  differences  are  significant .  From  a  cursory  look  at  the  available 

data,  the  performance  of  the  keCD  model  appears  to  provide  a  more  consistent 
behavior.  In  particular,  turbulence  stress  and  kinetic  energy  decrease  with  increasing 

convective  Mach  number  as  predicted  by  the  keCD  but  not  by  keCC. 

Relative  to  compressible  modeling,  a  recent  paper  by  Situ  and  Schetz  (6-91) 

considered  the  addition  of  the  fluctuating  term  -T  pV  to  the  shear  stress  where  a 
favorable  comparison  to  a  variety  of  shear  layer  and  boundary  layer  flows  was 
demonstrated.  Appendix  |  provides  details  of  this  turbulence  model  concept. 


3.4  Vortical/Curative  Modifications 

The  basic  k-e  model  was  developed  for  simple  shear  flows  (planar  shear  layers 
and  boundary  layers).  The  turbulence  field  is  assumed  to  be  isotropic  and  the  relating 
Boussinesq  eddy  viscosity  assumption  turbulent  stress  to  mean  strain  is  employed. 
For  flows  with  any  type  of  geometric  complexity,  corrections  are  required  to  enforce  the 
desired  behavior.  Corrections  are  required  for: 


(1)  axisymmetric  behavior 

(2)  stream  wise  curvature 

(3)  swirl 

(4)  streamwise  vorticity 

and  various  other  geometric  distortions  for  simple  planar  flow,  which  are  often  very 
specialized  and  problem  dependent.  The  corrections  range  from  those  being  entirely 
heuristic  to  those  attempting  to  be  generalized  and  thus  modeling  the  turbulent 
behavior.  All  the  latter  models  achieve  improved  behavior  via  modifications  to  the  e 
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equations  and  thus  they  directly  focus  on  the  length  scale  issue.  There  are  far  too 
many  corrections  in  usage  to  describe  in  this  report.  Rather,  several  will  be  briefly 
described  to  provide  an  understanding  of  how  such  corrections  are  made  and  how 

they  improve  the  performance  of  the  basic  k-6  model. 

The  first  corrections  to  the  k-e  model  were  made  in  the  early  ‘70’s  to  deal  with 
the  axisymmetric  problem,  namely,  the  inability  of  the  basic  model  to  analyze  both 

planar  and  axisymmetric  jets  with  the  same  constants.  To  have  the  k-e  turbulence 
model  agree  with  basic  low-speed,  axisymmetric  jet  data,  the  constants  must  be 
revised  or  made  to  vary  with  jet  parameters,  or,  additional  terms  must  be  added  to  the 
equations.  The  simplest  modification  entails  changing  C-j ,  from  a  value  1.43  for  planar 

to  1 .60  for  axisymmetric  jets.  This  forms  the  basic  for  some  of  the  earlier  work  with  k-e 
in  Great  Britain  in  analyzing  axisymmetric  rocket  plumes. 

The  next  group  of  corrections  make  the  coefficients  dependent  on  the  jet 
centerline  velocity  decay,  dUa/dx,  and  are  summarized  below. 

Launder,  et  al.  (1972)  ke1/ke2  models  provide  for 

C2  =  1.92-  aF,  a  =  0.067  (=0.053  for  ks2) 

=  .09  -  bF,  bsO.1 

where 


McGuirk  and  Rodi  (1977)  chose  the  coefficient  C-,  to  be 

Ci  =  1.14-  5.31  (42) 

Ucl  dx  '  ' 

while  Morse  (1977)  selected 

Ci  =  1.4  -  3.4  (43) 

In  the  above  axisymmetric  corrections,  U^,  is  the  jet  centerline  velocity  and  r1/2, 
is  the  distance  from  the  centerline  to  the  half-velocity  position. 

All  the  above  corrections  are  heuristic,  use  centerline  values  to  change  the 
turbulence  away  from  the  centerline,  and  are  not  extendible  to  more  generalized  flows. 
Pope  (3-78),  noting  these  issues,  introduced  a  more  generalized  extension  based  on 
vortex  stretching  arguments  which  modifies  the  dissipation  equation  by  adding  an 
extra  term,  namely: 
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where 


C3  X  e2/k 


(44) 


X  =  <°ij  “jk  Sy 


(45) 


with 


and 


Sij  =  i-k 
13  2  e 


3uj  9u,' 

— L  +  — I 


(3xj  9xjJ 


2  e  \9xj  3xjJ 
For  an  axisymmetric  jet  (with  no  swirl),  x  becomes: 

Y-irfcVr---Vi 

z  4  Ve/  Var  ax/  r 


(46) 


(47) 


(48) 


and  the  constant  C3  =  .79  was  found  to  match  jet  data  (axisymmetric  jet  into  still  air). 
This  vortex  stretching  correction  has  potential  generality  and  does  go  to  zero  for  planar 
flows.  However,  in  a  review  article  by  Launder  (1985)  the  universality  of  the  C3 

coefficient  is  questioned  based  on  calculations  of  coaxial  jets.  Hanjalic  and  Launder 
(1980)  suggest  that  a  correction  term  of  the  form: 


-  C3k 


9ui  din 
dxj  dxm 


Eijk  Sunk 


be  added  which  is  similar  to  Pope’s  in  appearance  but  acts  differently, 
different  and  this  term  is  also  operative  for  planar  shear  flows. 


(49) 

The  sign  is 


To  illustrate  the  importance  of  these  axisymmetric  correction  terms,  consider  a 
simple  low  speed,  axisymmetric  jet  with  a  10:1  velocity  ratio.  Figure  3  shows  the 

predicted  centerline  velocity  decay  and  jet  boundaries  obtained  using  the  basic  k-e 

model  and  the  k-e  model  with  the  Launder  CL  (centerline)  and  Pope  vortex  stretching 
corrections  (see  Dash,  6-91  and  Dash  et  al.,  9-91  for  details  and  discussions 
regarding  the  use  of  these  corrections  with  the  previous  compressibility  corrections 
provided  in  this  section).  As  exhibited,  inclusion  of  the  axisymmetric  correction  slows 
down  the  mixing  and  produces  a  significantly  larger  jet  size.  For  low  speed  jets  such 
as  this,  the  correction  produces  good  agreement  with  data. 


The  next  corrections  discussed  are  those  for  streamline  curvature.  While  both 
the  mean  flow  and  turbulence  model  equations  contain  a  number  of  curvature  terms 
arising  from  the  transformation  to  curvilinear,  surface-oriented  coordinates,  numerous 
investigators  have  demonstrated  that  additional,  curvature  corrections  terms  are 
required  to  account  for  the  strong  effects  of  curvature  on  turbulence  structure.  The 
analogy  drawn  by  Bradshaw  (1969)  between  curvature  and  buoyancy  has  been 
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utilized  by  most  investigators  as  the  basis  for  heuristic  corrections  to  algebraic  or  two- 
equation  turbulence  models.  Defining  the  curvature  parameters,  s  =  -KU/9U/9n),  a 
curvature  correction  to  the  ML  formulation  (Bradshaw,  1973)  is  given  by: 

Lp  (1  -  a  s)  (50) 

(1-s) 

where  Lq  is  the  planar  value  while  a  is  a  constant  (5  <  a  <  10).  The  treatment  has 

been  implemented  for  wall  jets  by  Foloyan  and  Whitelaw  (1976)  who  utilized  a 
complete  (inner/outer)  mixing  length  formulation.  For  the  present  near-wall  use  of  the 
mixing  length  formulation,  this  correction  will  only  be  required  in  situations  with  very 
large  curvature  since  the  near-wall  correction  to  Lq  is  negligible. 


Launder  and  coworkers  (3-77)  have  developed  a  curvature  correction  for  the  k-e 
model  which  utilized  a  single  empirical  coefficient,  Cc.  The  curvature  correction  is 
proportional  to  a  Richardson  number,  Ri,  based  on  the  turbulence  time  scale.  In  their 
formulation,  the  local  Richardson  number  is  given  by: 


Ri  = 


(51) 


and  the  C2  coefficient  of  the  e  equation  is  modified  as  follows: 


C2  =  1.92  (1  -  Q  Ri) 


(52) 


Values  of  Cc  of  about  0.2  have  yielded  optimal  predictions  for  a  variety  of  curved 
boundary  layer  flows  as  described  by  Launder  et  al.,  (3-77).  An  analogous  type  of 
curvature  correction  for  the  k-E  model  has  developed  by  Hah  and  Lakshminarayana 
(10-80).  They  have  modified  the  C-,  coefficient  of  the  e  equation  as  follows: 

Ci  =  1.43  (1  +CcRi)  (53) 

-  *  To  illustrate  the  marked  influence  of  these  streamwise  curvature  corrections,  we 
refer  to  the  work  of  Dash  et  al.,  (10-85  and  NASA  CP  2432,  1986)  for  wall  jets  on 
curved  surfaces.  The  case  simulated  the  experiment  of  Wilson  and  Goldstein  (9-76) 
(Res  =  13,200)  for  a  wall  jet  issuing  from  a  .615  cm  slot  flowing  over  a  circular  cylinder 
with  a  diameter  of  20.32  cm.  The  data  for  the  jet  half-radius  variation  is  compared  with 
a  series  of  predictions  in  Figure  4.  The  following  predictions  were  performed: 

(a)  Basic  k-e  (inner  Van  Driest) 
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(b)  k-e.  Launder  Correction 

As  above  the  k-e  curvature  correction  of 
Launder  using  Cc  =  .10  and  .20. 

(c)  As  above  with  k-e  curvature  correction,  eq.  (53), 
using  Cc  =  .12  and  .16,  Hah  and  Lakshminarayana 

(10-80) 

Comparing  the  data  for  this  case  with  that  for  the  planar  jet,  we  observe  a 
pronounced  effect  of  curvature  on  the  jet  growth  rate  (i.e.,  at  x/b  =  50,  the  planar  jet 

half-radius  is  about  4  times  larger).  Applications  of  the  basic  curvilinear  k-e  model 
grossly  underestimates  the  growth  rate  yielding  a  variation  comparable  to  that  of  the 

planar  case.  Launder  et  al.,  (3-77)  recommend  a  value  of  Cc  =  .20  for  their  k-e 
curvature  correction  but  state  that  Priddini  (1975)  had  found  a  value  of  Cc  =  .10  to 
work  best  for  wall  jets.  We  have  implemented  both  these  values  and  find  the  value  of 
.20  to  yield  best  agreement.  For  the  C,  curvature  correction  of  eq.  (53)  (using  the 

Launder  Ri  definition),  the  value  of  Cc  =  .16  is  seen  to  yield  good  agreement  with  the 
data  while  the  value  of  Cc  =  .12  underestimated  the  measured  growth  rate. 

Comparable  conclusions  can  be  drawn  from  the  comparison  of  predictions  for 
maximum  velocity  decay  also  depicted  in  Figure  4.  The  predictions  with  the  Launder 
correction  with  Cc  =  .20  and  the  Cj  correction  with  Cc  =  .16  are  quite  comparable  and 

fit  within  the  data  band  based  on  the  experiments  of  Wilson  and  Goldstein,  and  other 
investigators.  The  other  k-e  predictions  underestimate  the  velocity  decay  rate. 

Predicted  profiles  of  the  streamwise  velocity  at  0  =  90°  (with  Cc  =  .20,  Launder 
and  Cc  =  .16  Eq.  (53)  are  indistinguishable  and  agree  favorable  with  the  data  of 

Wilson  and  Goldstein  which  is  invariant  in  similarity  coordinates  for  0  =  45°.  The  case 
itself  is,  however,  nonsimilar  as  seen  by  the  comparisons  with  turbulent  shear  stress 
data  depicted  in  Figure  5.  The  data  exhibits  a  significantly  greater  peak  shear  stress 
level  at  1 80°  than  at  90°  as  do  the  predictions.  The  predictions  at  90°  <Cc  =  .20, 

Launder  and  Cc  =  .16  with  eq.  (53)  are  in  close  agreement  with  each  other  and 

somewhat  overestimate  the  shear  stress.  Those  at  180°  differ  substantially  with  the 
Launder  model  overestimating  and  the  Hah  and  Lakshminarayana  model 
underestimating  peak  shear  stress  levels.  Figure  5  also  compares  the  predicted 
streamwise  variation  of  peak  turbulent  shear  stress  with  data.  The  two  predictions  are 
quite  comparable  up  to  70°  with  the  eq.  (53)  model  yielding  a  decreased  rate  beyond 
70°  while  the  Launder  model  continues  in  an  essentially  linear  fashion.  Discounting 
the  data  point  at  1 80°,  the  trend  provided  by  eq.  (53)  appears  to  be  more  consistent; 
moreover,  the  model  yields  values  of  peak  turbulent  shear  stress  which  agree  quite 
closely  with  the  data  of  0  <  0  <  150°. 
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4.0  MULTI-ZONE  NAVIER-STOKES  CODE 


The  code  selected  for  turbulence  model  assessment  studies  is  the  PARCH 
Navier-Stokes  (NS)  code.  PARCH  is  an  outgrowth  of  the  Ames  Research  Code  (ARC) 
of  Pulliam  (3-84  and  1-85)  and  the  propulsive  extension,  PARC,  that  Cooper  (6-90) 
developed  at  AEDC.  Like  the  AFWAL  PNS  code,  PARCH  embodies  over  ten  years  of 
continuous  development,  upgrade,  and  support  by  the  government. 
ARC/PARC/PARCH  represent  the  most  popular  and  widely  used  NS  codes  in  the 
United  States.  The  PARCH  version  has  unified  the  features  of  ARC  and  PARC  and 
includes  unique  “patched  grid”  capabilities  developed  at  AEDC  to  treat  complex 
geometries  using  simple  grldding  procedures.  With  patched  gridding,  boundary 
conditions  can  be  applied  along  any  mapped  grid  lines  of  the  computational  domain, 
not  just  the  bounding  domain.  This  facilitates  dealing  with  fins,  struts,  cavities,  etc. 

PARCH  is  an  implicit  code  which  is  requisite  for  the  efficient  treatment  of 
turbulent  flows  near  walls  where  grid  spacing  is  small  and  would  severely  limit  time- 
steps  if  explicit  methods  were  utilized.  Separate  2D/AXISYMMETRIC  and  3D  versions 
are  operational  which  have  been  specialized  to  various  flow  problems  shown  in  Table 
VII.  Versions  of  PARCH  are  operational  with  finite-rate  chemistry  and  multi-phase  flow. 
Multi-zone  versions  are  available  where  each  zone  has  internal  patching  capabilities 
as  well.  PARCH  has  been  applied  to  aerodynamic  flows  (hypersonic  vehicles,  tactical 
missiles,  circulation-control  airfoils),  to  propulsive  flows  (rocket  nozzles,  gas  turbines, 
scramjet  inlets/combustions/nozzles),  jets/plumes/wakes  (aircraft  jets,  rocket  plumes, 
RV  wakes),  and  laser  flows.  Moreover,  it  is  presently  being  used  as  a  research  tool 
under  NASA/LaRC  support  to  study  advanced  turbulence  models  for  jets  (in  support  of 
the  high  speed  civilian  transport  program). 

Turbulence  models  have  been  incorporated  into  PARCH  in  both  a  loosely- 
coupled  manner  and  a  strongly-coupled  fashion.  The  PARCH/Gas  Turbine  version 
presently  contains  the  following  turbulence  models: 

(1 )  k-e  model,  high  Re  version 

(2)  compressibility-corrected  variants  of  the  above  with  the  corrections 
of  Dash  (8-75),  Sarkar  et.  al  (12-89),  and  Zeman  (2-90) 

(3)  near-wall  extensions  including  van-Driest  two-layer  coupling,  the 
Chien  (1982)  low  Re  version,  and  the  Rodi  (1-91)  one-equation 
hybrid  approach 

(4)  geometric  corrections  for  axisymmetric  effects  (Launder,  1985  and 
Pope  (3-78)  corrections),  for  swirl  and  streamwise  curvature 
corrections  (Penn  State/Lakshminarayana,  12-86  and  Imperial 
College) 

(5)  k-w  model,  high  Re  yersion 
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Table  VII  (cont.).  SPECIALIZED  VERSIONS  OF  PARCH  CODE 
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PARCH  has  been  applied  to  a  very  large  data  base  of  fundamental  flows  as  well 
as  to  many  complex,  three-dimensional  flows  as  demonstrated  by  Dash  et.  al  (6-91,  1- 
91,  and  10-89),  York  et.  al  (6-89  and  11-89),  and  Sinha  et.  al  (7-88  and  7-89). 

The  selection  of  PARCH  as  a  generalized  tool  to  study  the  behavior  of 
turbulence  medels  in  high  speed  flows  is  motivated  by  many  new  practical  factors. 
With  current  generation  missiles  and  spacecraft  having  highly  integrated  propulsive 
systems,  one  cannot  separately  study  aerodynamic  and  propulsive  flows. 

The  exhaust  flow  from  integrated  propulsive  systems  can  wet  large  protions  of 
the  aerodynamic  surfaces.  This  flow  can  be  chemically-reacting  and  for  missiles  with 
solid-propellant  rockets,  will  be  multi-phase.  With  lateral  control  jets  replacing 
conventional  aerodynamic  control  surfaces  on  many  missiles,  the  ability  to  analyze 
such  flowfields  is  felt  to  be  requisite.  PARCH  which  has  chemical  and  multi-phase 
capabilities,  is  the  only  3D  Navier-Stokes  Code  now  being  used  in  a  design 
environment  to  support  the  development  of  varied-missile  systems  with  integrated 
propulsive  systems  and  missiles  with  lateral  control  jets. 

Another  selection  criterion  for  PARCH  is  the  validation  of  the  code  with  its 
included  turbulence  models.  This  has  occured  over  a  five  year  period  starting  with  the 

inclusion  of  the  basic  high  Re  version  of  k-e  and  checkout  for  varied  high  speed  free 

shear  flows,  and  versions  of  k-e  coupled  to  various  near-wall  models.  It  takes  many 
man-years  of  effort  to  incorporate  an  advanced  turbulence  model  into  a  Navier-Stokes 
code  and  get  it  to  work  properly  in  a  broad  spectrum  of  flow  problems.  No  other 
Navier-Stokes  code  has  the  experience  base  of  PARCH  with  advanced  turbulence 
models. 

Finally  PARCH  is  now  used  as  a  “standard”  tool  by  many  different  communities 
dealing  with  aerodynamic,  propulsive,  signature  and  acoustic  problems.  PARCH  is  a 
standard  component  of  the  SDIO  system  of  codes  for  missile  detection,  the  SPIRITS 
system  of  codes  for  aircraft/helicopter  detection,  the  HSCT  system  of  codes  (MICOM) 
NSWC,  etc)  for  missile  design,  etc.,  etc.  All  these  communities  are  dealing  with  the 
turbulence  problem  in  a  practical  manner  and  their  experience  with  PARCH  will 
provide  significant  contributed  value  to  envisioned  work  by  Wright-Labs  with  this  code. 
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5.0  TURBULENCE  MODELING  ROADMAP 


In  order  to  accommodate  a  wide  range  of  practical  engineering  problems  using 
a  generic  turbulence  closure  model,  a  “building  block”  approach  will  be  required.  This 
approach  consists  of  evaluating  a  variety  of  turbulent  boundary  layer  flows  relative  to 
corresponding  validation  cases  which  can  be  used  to  calibrate  closure  constants  (or 
functional  coefficients)  together  with  appropriate  damping  terms.  The  methodology 
would  commence  with  simple  unit  level  problems  and  progress  to  more  complex  flows 
with  combined  effort  as  the  experience  based  is  enhanced.  Table  VIII  shows  the 
descending  order  of  flow  status  that  should  be  considered.  The  end  product  of  the 
methodology  will  provide  a  matrix  (table)  of  closure  coefficients  and  damping  functions 
that  can  be  used  in  the  proposed  hybrid  turbulence  model  that  will  allow  the  user  to 
“dial-a-model”  as  required  to  solve  the  problem  at  hand.  Several  techniques  that 
can  be  used  to  develop  the  building  block  methodology  will  be  addressed  in  this 
section  as  well  as  the  identification  of  a  fundamental  database  for  validation. 

The  building-block  matrix  shown  in  Table  VIII  will  consist  of  compressible 
boundary  layer  and  jet  flowfields  that  are  2D  or  axisymmetric.  While  several 
successes  have  been  demonstrated  for  unit  problems,  coupling  of  unit  level  problems 
with  compressibility  has  not  been  established.  For  a  generic  turbulence  model  to 
analyze  fundamental  wall  bounded  and  jet  flowfields,  the  turbulence  modeling  in  each 
flow  state  must  be  validated  in  a  reliable  manner.  Consequently,  a  systematic 
development  from  the  simplest  flow  to  the  more  complex  flows  will  be  required.  The 
emphasis  of  our  approach  will  be  on  boundary  layer  type  flows.  A  similar  approach  for 
jet  flowfields  has  been  proposed  to  NASA  (Dash  et  al,  10-91)  that  would  complement 
proposed  work  for  a  Phase  II  program. 


5.1  Wall  Roughness 

Appendix  A  addressed  issues  of  wall  roughness  in  boundary  layer  flow.  The 
modification  to  the  C-S  model  to  incorporate  the  van-Driest  local  vortex  generating 
factor  was  shown  to  be  inadequate  since  the  factor  did  not  allow  for  roughness 
elements  greater  than  the  sublayer  thickness.  Since  the  near  wall  velocity  distribution 
can  be  dependent  upon  many  geometric  variables  that  define  surface  texture,  the 
concept  of  Rotta  (1962)  will  be  adapted.  Here,  the  influence  of  surface  roughness  is 
felt  only  close  to  the  wall  such  that  the  law  of  the  wall  relations  could  be  used  by 

shifting  the  reference  plane  by  some  distance  (-A  y+)  with  respect  to  mean  roughness 
level.  To  accommodate  the  no-slip  condition  at  the  new  surface,  the  reference  plane 
would  be  required  to  move  in  the  direction  opposite  to  the  mean  flow  at  a  rate  that 
depends  upon  the  shift. 

Cebeci  and  Chang  (7-78)  used  the  Rotta  criteria  in  an  algebraic  mixing  length 

model  but  allowed  for  the  no-slip  condition  to  be  maintained  at  y+  =  0.  Chang  (5-71) 
used  reference  plane  slip  velocity  but  did  not  allow  for  a  shift  in  the  reference  plane 
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Table  viii.  Building-Block  Turbulence  Modeling 
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Wave/Shear  .  Unbalanced  Press. 

Layer  Interaction  Above  &  Wave/Shear 

Layer  Interactions 


position.  Laganelli  (1-75)  as  well  as  Krogstad  (6-91)  considered  both  aspects  of  the 
Rotta  concept.  The  latter  paper  provided  a  thorough  examination  of  the  damping 
requirements  and  showed  excellent  comparisons  to  incompressible  data.  In  both 
cases,  the  authors  chose  to  characterize  the  reference  surface  by  a  single  roughness 
length  scale,  namely,  the  equivalent  sandgrain  roughness  height  which  is  easily 
obtained.  In  a  recent  investigation,  Laganelli  and  Scaggs  (6-90)  reviewed  the  single 
sandgrain  roughness  length  scale  criteria  of  Dirling  (6-73),  which  was  developed  from 
incompressible  data,  and  demonstrated  that  the  criteria  could  be  used  for 
compressible  flow  conditions.  The  authors  (Laganelli  and  Scaggs,  1-86)  determined 
that  a  unit  Reynolds  number  (velocity  scale)  was  required  in  the  equivalent  sandgrain 
roughness  parameter.  The  authors  also  reviewed  the  reference  plane  shifting  criteria 
using  a  boundary  layer  code  and  recommended  a  technique  to  define  the  reference 
plane  position.  A  comparison  to  compressible  data  showed  good  agreement. 

When  the  roughness  elements  just  protrude  through  the  viscous  sublayer, 
turbulence  production  can  be  generated  by  the  individual  roughness  elements  where 
the  eddy  length  is  on  the  order  of  the  element  size.  In  regions  away  from  the  wall,  the 
impact  of  the  vortex  generating  roughness  elements  could  have  an  effect  on  local 
equilibrium  (i.e.  production  is  approximately  equal  to  dissipation)..  However,  with  the 
use  of  two-equation  models,  the  field  equations  should  provide  insight  into  near-wall 
modeling  requirements  especially  since  the  kinetic  energy  is  carried  to  the  wall. 

Finally,  it  should  be  noted  that  Arora  et.  al  (11-82)  used  the  technique  of  Cebeci 

and  Chang  in  a  two-equation  k-e  model  (using  the  wall  function  method)  and  the 
assumption  of  equilibrium  away  from  the  wall.  Comparison  of  the  method  to 
incompressible  and  compressible  data  showed  good  agreement  which  could  render 

the  no-slip  condition  y+  =  o  as  opposed  to  the  reference  plane  as  not  important. 

The  velocity  distribution  in  the  wall  region  given  by  Krogstad  (6-91)  is 


u 


2&y*  _ 

Jo  H^l*<2rDR)* 


(54) 


where  i+  is  the  dimensionless  mixing  length  and  DR  is  a  roughness  damping  function 
(given  in  the  bracket  of  Eq.  A-12).  Since  DR  is  limited  to  DR  <  1  (as  given  in  Eq.  A-12), 
Krogstad  proposed  a  shift  such  that 
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An  examination  of  the  above  shows  the  damping  term  to  be  similar  to  that  given  by 
Eqs.  (A-13)  and  (A-14).  However,  Eq.  (55)  appears  to  have  correct  asymptotic 
behavior  in  both  the  linear  region  of  the  sublayer  as  well  as  the  log-law  (slope)  region 
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that  correctly  accounts  for  the  velocity  shift.  The  form  given  by  Eqs  (A-13)  and  (A-14) 
used  an  average  value  in  the  sublayer  region  and  does  not  have  the  asymptotic 
characteristics  of  Eq.  (55). 

This  work  should  provide  excellent  near  wall  modeling  techniques  to 
incorporate  into  the  proposed  hybrid  k-e  turbulence  model. 


5.2  Mass  Transfer 

Appendix  A  also  describes  the  basic  methodology  that  has  been  used  by  the 
modelers  to  treat  mass  transfer  effects.  Essentially,  corrections  have  been  made  to  the 

damping  coefficient  such  that  A+  =  A+  (vw).  These  corrections  have  been  adapted  into 

algebraic  eddy  viscosity  models  as  well  as  two-equation  models.  For  the  k-e  type, 
modifications  are  usually  made  to  damping  coefficients  in  the  expontial  term  in  wall 
functions,  van  Driest  bridging  techniques,  and  low  Reynolds  number  near-wall 

techniques.  The  k-cotype  have  also  featured  mass  transfer  modeling  using  law  of  wall 
relations  coupled  to  boundary  conditions  using  a  perturbation  solution  in  the  near  wall 
region.  However,  all  of  the  above  methods  have  shown  success  for  incompressible 
flow  conditions  only.  To  the  authors  best  knowledge,  no-one  has  demonstrated  the 
ability  to  predict  mass  transfer  effects  subject  to  compressible  flow. 

Rajendran  and  Laganelli  (12-89)  considered  the  arguments  of  van-Driest  that 
the  presence  of  the  wall  can  modify  the  universal  constant  (Karman  constant)  such  that 

k  =  k  [1-exp  (-y/A)] 

or  that  the  mixing  length  must  be  charged  to 

i  =  xy  [1-exp  (-y/A)] 

An  examination  of  compressible  mass  transfer  data  of  Martellucci  et  al  (4-74) 
and  Danberg  (1-67)  shows  that  the  slope  in  law  of  the  wall  coordinates  deviates  from 
the  non-blowing  value  and  changes  with  increasing  blowing.  Accordingly,  Rajendran 
and  Laganelli  suggest  that  the  damping  term  he  modified  in  the  same  spirit  as  wall 
roughness  in  as  much  as  wall  blowing  has  the  opposite  effect  as  roughness  in  law  of 
the  wall  coordinates.  For  this  condition,  the  damping  can  be  expressed  as 

[l -exp(-y*/A*)\  +  ${y*/A*)a  ex p  [Bu/BUer]  (56) 


where 


P  =  p(X)  =  (pv)w/(pu)e 


(57) 


49 


Bu  =  (X/St) 


(58) 


=  exp  [2.7726  (t/-0.5)  ] 


W  -  (Tw/Te)~‘12S  +  0.125W, 


(59) 


In  the  above,  the  critical  blow-away  parameter  Bucr  was  determined  by 
Laganelli  et.  al.  (4-75)  to  vary  between  a  value  of  4  (incompressible)  to  64 
(compressible)  using  data  from  a  several  facilities  over  a  range  of  Mach  number  and 
wall  temperature  ratio.  The  technique  was  used  in  the  Baldwin-Lomax  algebraic  eddy 

viscosity  model  and  the  constants  (3  and  n  were  determined  from  the  data  of 
Martellucci  et  al  (4-74).  This  initial  result  showed  promise  and  will  allow  for  a 
technique  to  be  further  evaluated  in  a  Phase  II  program. 

With  the  ability  to  provide  more  physics  to  the  damping  coefficient  A+  using  the 
kinetic  energy  equation,  together  with  the  heuristic  arguments  developed  for  near  wall 
functions,  we  believe  that  the  ability  to  characterize  mass  transfer  can  be 
accomplished.  This  can  be  demonstrated  by  comparing  the  exponent  terms  from  the 
damping  functions  of  van  Driest  and  Norris-Reynolds  i.e. 
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which  leads  to  the  definition  of  the  damping  coefficient  to  be  A+  =  Ux  /C3  k1/2.  Here, 

one  has  the  added  advantage  of  the  kinetic  energy  which  is  provided  by  the  field 

equation  and  should  provide  better  physics  for  determining  characteristics  of  A+  =  A+ 
(vw)  as  opposed  to  data  correlations. 

5.3  Compressibility 

Section  3.3  addressed  issues  concerning  compressible  modifications  to  the 

standard  k-e  incompressible  field  equations.  Specifically,  investigations  were  made 
by  Dash  et  al  (6-91  and  9-91),  Viegas-Rubesin  (6-91),  and  Wilcox  (6-91)  that 
considered  adaption  of  the  Sarkar  et  al  (12-89  and  6-90)  and  Zeman  (2-90) 
compressible-dissipation  models  to  the  k-e  and  k-co  two-equation  turbulence  models, 
respectively.  While  Dash  et  al  and  Viegas-Rubesin  focussed  on  mixing  layers,  Wilcox 
considered  both  mixing  layers  and  wall  bounded  flows.  The  authors  noted  some 
variation  in  the  compressible  models. 
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Several  observations  can  be  drawn  from  the  work  of  the  above  investigators 
which,  when  supplemented  by  the  work  of  Horstman  (6-91 )  provide  a  practical  path 
that  can  be  followed  to  enhance  the  prediction  of  compressible  shear  flows.  Since 
section  3.3  focussed  on  free  mixing  layers,  the  immediate  near  path  should  consider 
the  impact  of  these  compressible  modifications  to  the  turbulence  two-equation  models 
in  analyzing  wall  bounded  flows. 

When  the  Sarkar  et  al  and  Zeman  compressible  terms  were  added  to  the  k-e 
and  k-co  models,  a  comparison  with  free  shear  layer  data  showed  that  the 
compressibility  effects  on  spreading  rate  was  reasonably  predicted.  However  when 
applied  to  adiabatic,  flat  plate  compressible  boundary  layer  flow  data  (M  <  5),  the 
compressible  corrections  tended  to  underpredict  skin-friction  data  while  the 
unmodified  turbulence  two-equation  model  predicted  this  data  quite  well.  Wilcox 
modified  the  coefficients  in  the  Sarkar  model  to  force  agreement  with  the  skin-friction 
data  and  subsequently  used  the  modified  version  to  analyze  the  free  shear  data  base. 
Wilcox’s  experiences  indicated  that  the  compressibility  modifications  can  not  be 
arbitrarily  extended  to  the  near  wall  region,  even  for  simple  parallel  shear  flows.  Initial 
application  of  the  Sarkar  compressible  modifications  to  separated  flow  conditions 
(M<3)  provided  a  larger  separation  length  than  the  unmodified  models  with 
significantly  better  agreement  with  data  (Dash,  9-91). 

Wilcox  (6-91),  demonstrated  that  the  Sarkar  and  Zeman  models  were 
inaccurate  in  the  near  wall  region  and  developed  an  improved  model  that  combines 
the  Zeman  extended  lag  effect  with  the  Sarkar  dilatational  dissipation  effects  on  the 
turbulent  Mach  number.  When  the  hybrid  model  was  applied  to  the  limited  data  base 
of  mixing  layers  (Mc<2),  adiabatic  flat-plate  boundary  layers  (M<5),  and  separated 

flows  (M<3),  data  trends  were  correctly  predicted.  The  k-co  equations  were  used  for 
this  examination. 

Horstman  (6-91)  used  the  standard  k-e  model  and  two  modifications  of  this 
model  and  compared  the  results  to  a  range  of  supersonic/hypersonic  separated  flows. 

The  modified  k-e  models  consisted  of  the  standard  k-e  equations  and  a  one-equation 

model  (hybrid,  Rodi,  1-91)  for  length  scale  in  the  near  wall  region  and  the  standard  k-e 
model  with  additional  compressible  terms  (as  provided  from  Farve  averaging  by 

Rubesin,  6-90).  The  standard  k-e  tended  to  underpredict  the  separation  length  and 
overpredict  skin-friction  and  heat  transfer.  The  Rodi  hybrid  model  tended  to  increase 
the  separation  zone  while  underpredicting  skin-friction  and  heat  transfer  in  the 

reattachment  region.  The  k-e  model  with  additional  compressible  terms  tended  to 
provide  the  same  trends  as  the  hybrid  model  but  was  difficult  to  implement  readily.  An 
examination  was  made  to  determine  the  impact  of  individual  effects  of  compressible 
modifications  and  length  scale  (Rodi)  modifications.  It  was  shown  that  both 
compressible  and  length  scale  modifications  increased  the  separation  zone  length 
while  length  scale  alone  changed  skin-friction  and  heat  transfer  in  the  reattachment 
region. 
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Based  on  such  assessments,  our  building  block  approach  will  consist  of 
“carefully”  combining  improvements  for  free  shear  layer  flows  and  boundary  layer 

flows  using  the  standard  k-e  model  with  compressible  modifications  of  Dash  (9-91)  as 
the  baseline  model.  Near-wall  behavior  will  be  explored  using  low  Reynolds  number 
modifications  (with  compressible  correction)  and  the  Rodi  hybrid  length  scale 
modifications  (also  with  appropriate  compressible  modifications). 

5.4  Pressure  Gradient  (Curvature) 

Section  3.4  addressed  issues  concerning  vortical-curvature  modifications. 
These  modifications  include  corrections  for  axisymmetric  behavior,  streamwise 
curvature,  swirl,  and  streamwise  vorticity.  As  noted,  these  modifications  have  been 

made  to  the  constants  in  the  e  equation  of  the  k-e  two-equation  model. 

Consistent  with  our  building  block  approach,  a  systemetric  evaluation  of  varied 
geometric  corrections  will  be  required.  It  is  imperative  to  couple  the  assessment  of  the 
varied  geometric  corrections  with  the  assessment  of  the  compressibility-correction 
(both  influence  the  length  scale),  and,  with  the  near-wall  formulation.  The  geometric 
corrections  are  generally  formulated  for  low  speed  flows  and  they  may  not  work 
without  recalibration  for  higher  speed  flows  even  with  the  compressibility-correction 
implemented. 

5.5  Compressible  Flow  Database 

The  previous  sub-sections  considered  issues  and  methodologies  that  have 
been  used  to  predict  attached/separated  boundary  layer  and  jet  flow  fields.  Moreover, 
several  of  the  methodologies  have  been  shown  to  predict  favorable  and  adverse 
pressure  gradient  flows  subject  to  compressible  flow  conditions.  Our  objective  has 
been  the  development  of  a  generic  turbulence  closure  model  capable  of  providing 
near  term  solutions  to  a  wide  range  of  engineering  problems.  Accordingly,  we  have 

identified  the  k-e  two-equation  model  with  variations  to  meet  these  near  term  goals. 

Inparticular,  this  will  involve  the  hybrid  concept  that  consists  of  the  standard  k-e  model 
and  a  one-equation  model  for  length  scale.  Variations  will  require  a  systematic 
investigation  to  provide  compressible  corrections,  curvature  corrections,  roughness 
and  blowing  corrections. 

A  database  will  be  required  to  provide  this  systematic  investigation.  The 
Stanford  conferences,  NASP  symposiums,  and  recent  AIAA  meetings/workshops  have 
established  a  formidable  list  of  test  results  which  can  be  used  to  develop  validation 
cases.  Table  IX  lists  a  database  of  compressible  flow  experiments  for  boundary  layer 
and  jet  flowfields.  This  preliminary  database  is  consistent  with  the  building  block 
approach  described  in  Table  VIII.  While  several  of  these  references  provide  data  for 
unit  flow  problems  (mass  transfer,  for  example),  other  define  experiments  for  several 
unit  level  problems.  The  database  allows  development  for  simple  flat  plate,  adiabatic 
parallel  shear  flows  to  complex  shock/boundary  layer  interactions  and  subsequent 
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Table  IX.  SOURCES  OF  COMPRESSIBLE  FLOW  DATABASE 
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lources:  STANFORD  II  (Kline  et  al.,  1982) 

Martellucci  et  al.  (PTP-S,  3-81) 
Bushnell  (1-91) 


coupling  of  flow  states.  Most  of  these  test  cases  (sans  mass  transfer  and  rough  wall 
data)  have  been  used  by  previous  modelers  using  conventional  techniques  and  in 
some  cases  with  modifications  to  standard  techniques.  The  quality  of  this  database 
has  essentially  been  judged  and  its  limitations  defined. 

Finally,  Table  X  shows  how  the  building  block  approach  can  be  used  to  develop 
a  matrix  for  a  turbulence  “dial-a-model”  to  suit  specific  classes  of  engineering 
problems.  It  is  our  position  to  use  the  hybrid  versions  of  the  k-e  turbulence  model  with 
appropriate  compressible  corrections  as  the  baseline  model.  Essentially  this  table 
would  be  structured  having  established  the  ability  to  predict  compressible  zero- 
pressure  gradient,  adiabatic  and  cold  wail  conditions  for  two- 
dimensional/axisymmetric  flow  fields. 

As  shown  in  Table  X  and  described  in  the  text,  near-wall  length  scales  will 
require  modifications  to  include  mass  transfer,  roughness,  and  curvature.  These 
modifications  will  require  adjustments  to  constants  (which  can  become  functional 
coefficients)  and/or  adjustments  to  a  reference  plane  which  requires  the  no-slip 

conditions  at  this  plane.  As  noted  by  Rotta,  this  requires  a  shift  (-A  y+)  relative  to  near 

roughness  level  and  +A  y+  (Rajendran  and  Laganelli)  for  mass  transfer.  For  curvature 

these  adjustments  can  be  made  to  either  length  scales  or  the  constants  in  the  k-e 
transport  equations. 

A  potential  variation  to  the  above  would  be  to  provide  a  functional  form  for  the 
constants,  C-j ,  C2,  and  o£  to  include  wall  effects.  We  have  seen  that  low  Reynolds 

number  forms  to  the  k-e  equations  modify  the  constants  C-j  and  C2  to  include  damping 
terms  (mass  transfer  and  roughness)  or  velocity  gradients  of  mean  flow  for  curvature. 

We  believe  that  the  baseline  compressible  model  together  with  the  extensive 

experience  base  of  k-e  modelers  to  predict  a  wide  range  of 
incompressible/compressible  flow  engineering  problems,  will  provide  the  insight  to 
complete  the  matrix  of  Table  X.  However,  this  will  require  the  systematic  building  block 
approach  and  database  described  above. 


Table  X.  Building-Block:  Dial-A-Model  Matrix  for  Turbulence  Closure  Near-Wall  Requirements 


Parameters 


6.0  ASSESSMENT  OF  PHASE  I  PROGRAM  AND  RECOMMENDATIONS 

FOR  PHASE  II  EFFORT 

6.1  Assessment  of  Phase  I  Program 

A  through  review  has  been  made  to  idenitfy  a  turbulence  closure  model 
capable  of  treating  a  wide  range  of  practical  engineering  problems.  The  review  has 
focussed  on  high  Mach  number  compressibility  effects  for  simple  turbulent  flows- 
boundary  layers  as  well  as  free  hear  layers.  Consideration  was  also  given  to  wall 
effects  such  as  roughness,  blowing,  and  curvarture  as  well  as  flows  experiencing 
separation  and  interactions.  The  process  of  selecting  a  turbulence  model  also 
considered  adaptability  to  a  multizonial  Navier-Stokes  solver  with  a  long  history  in  the 
industrial  community. 


It  was  determined  that  a  hybrid  k-e  model  had  the  potential  to  meet  the  above 
objections  for  near  term  resolution  to  engineering  problems.  Specifically,  the  hybrid 

model  consists  of  the  standard  two-equation  high  Reynolds  number  k-e  turbulence 
model  for  treating  the  outer  layer  and  a  one-equation  model  (with  correct  length 

prescription  and  the  e  of  the  k-e  two-equation  model)  to  treat  the  inner  viscous  layer. 
The  hybrid  model  which  was  introduced  by  Rodi  for  incompressible  flows  has  been 
successfully  applied  to  analyze  compressible  flows  including  complex  shock/boundary 
layer  interactions.  The  hybrid  model  was  chosen  for  a  number  of  reasons,  namely: 

•  k-e  two-equation  model  is  the  most  widely  used  turbulence  model 
in  the  scientific  community 

•  it  has  an  extensive  experience  base  to  treat  a  wide  range  of 
engineering  problems 

•  use  of  the  one-equation  model  to  treat  the  length  scale  in  the  near¬ 
wall  region  eliminates  numerical  problems  associated  with  the  e 
equation  in  the  k-e  two-equation  model 

•  one-equation  models  have  demonstrated  success  to  provide 
engineering  solutions  for  flows  subject  to  serve  adverse  pressure 
gardients  and  has  demonstrated  success  to  provide  engineering 
solutions  for  hypersonic  flows  subject  to  complex  shock/boundary 
layer  interactions 

•  provides  details  of  the  turbulence  structure  as  well  as  wall 
conditions 

•  k-e  transport  equations  have  demonstrated  success  in  adapability 
to  industrial  codes 
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Several  turbulence  models  were  reviewed  that  included  zero-equation  type 
where  both  length  and  velocity  scales  are  algebraically  prescribed,  one-equation 
models  where  the  length  scale  is  prescribed  and  the  velocity  scale  is  determined  by  a 
kinetic  energy  transport  equation,  and  the  two-equation  models  where  both  length  and 

velocity  scales  are  determined  by  field  equations;  the  hybrid  k-e  being  in  the  latter 
class.  Variation  in  these  models  that  have  shown  success  include  the  one-half 
equation  model  of  Johnson  and  King,  the  one-equation  models  of  Goldberg  and 
Baldwin  and  Barth,  and  a  number  of  variations  to  the  standard  high  Reynolds  number 

k-e  model.  Recent  success  has  also  been  reported  by  Situ  andSchetz  who  included 
extra  terms  in  the  basic  shear  stress  equation  to  treat  compressibility  effects. 

Finally,  Dash  et  al,  Wilcox,  Horstman,  and  Viegas  and-Rubesin  considered  the 
impact  of  compressibility  with  modifications  to  the  two-equations  turbulence  models. 
Horstman  considered  the  extra  compressible  terms  suggested  by  Rubesin  (cascadsed 

from  Favre  averaging)  and  compared  this  result  to  the  standard  k-e  and  the  hybrid  k-e 
noted  above.  It  was  determined  that  both  the  hybrid  k-e  and  the  extra  compressible 

terms  k-e  model  provided  better  resolution  to  interacting  flows  where  the  former  gave 
better  control  on  separation  zone  length  and  magnitude  of  wall  heat  transfer.  The  work 
by  Dash,  Wilcox,  and  Viegas-Rubesin  considered  adaption  of  the  Sarkar  et  al  and 
Zeman  compressible  terms  to  the  two-equation  models.  It  was  determined  that  the 
compressible  terms  provided  the  reduction  in  the  spread  rate  of  mixing  layers;  but 
when  applied  to  boundary  layer  flows  (Wilcox)  underpredicted  results  compared  to 
unmodified  turbulence  models.  However,  the  Wilcox  model  is  not  hybrid  and  our  use 
of  k-e  with  the  Rodi  near  wall  formulation  may  remedy  this  deficiency. 

Recognizing  the  potential  of  the  hybrid  k-e  turbulence  model  and  the  issues 
concerning  compressibility,  a  building  block  approach  has  been  identified  that  can 

provide  a  matrix  of  turbulence  closure  parameters  that  can  be  used  with  the  hybrid  k-e 
model  to  solve  a  wide  range  of  practical  engineering  problems.  This  matrix  will 
provide  the  user  with  a  “dial-a-model”  menu  capable  of  analyzing  boundary  layer  and 
free  shear  layer  flows.  In  order  to  develop  this  matrix  through  a  building  block 
approach,  a  database  and  a  multizone  Navier-Stokes  solver  were  identified.  The 
database  consists  of  attached/separated  boundary  layer  flows  as  well  as  plane  free 
shear  axisymetric  free  jet  flows.  This  compressible  flow  database  has  been  used  by 
the  scientific  community  and  has  been  judged  for  quality  and  limitations  for  validating 
models. 

The  multizone  Navier-Stokes  solver  identified  for  the  building-block 
investigation  is  the  PARCH  code  which  is  an  outgrowth  of  the  Ames  Research  Code 
(ARC)  and  the  propulsive  extension,  PARC,  developed  at  AEDC.  PARCH  embodies 
over  ten  years  of  continuous  development,  upgrade,  and  support  by  various  DoD 
agencies.  ARC/PARC/PARCH  represent  the  most  popular  and  widely  used  NS  codes 
in  the  United  States.  A  simplified  version  of  the  PARCH  code  that  does  not  include 
chemistry  is  recommended  for  the  development  of  the  turbulence  closure  matrix. 
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6.2  Recommendations  for  Phase  II  Effort 


A  Phase  II  program  is  recommended  that  will  focus  on  boundary  layer  type 
flows.  Similar  efforts  are  currently  being  proposed  to  NASA  (by  Sciences  Application 
International  Corporateion)  to  investigate  free  shear  and  axisymmetric  jet  flows.  Both 
of  these  programs  would  compliment  each  other  and  could  be  developed  in  a  con¬ 
current  fashion  for  technical  interchange  and  cost  effectiveness. 

Prior  to  excersing  the  building  block  methodology,  it  will  be  necessary  to  resolve 
the  compressibility  issues  identified  in  the  Phase  I  investigation.  The  database  will  be 
used  together  the  concepts  proposed  in  this  investigation;  namely  to  include  the 
Sarkar  and  Zeman  compressible-  dilatation  terms,  or  some  modification  of  the  two, 

into  the  hybrid  k-e  model.  It  is  also  recommended  that  a  low  Reynolds  number  version 

of  the  k-e  model,  such  as  Chien,  be  modified  in  the  same  fashion  for  comparison. 
Having  established  a  baseline  compressible  turbulence  model,  the  building  block 
methodology  can  proceed  from  simple  to  more  complex  flows  to  treat 
attached/separated  flows  including  roughness,  blowing,  and  curvature.  Unit  level 
problems  will  be  investigated  to  develop  the  “dial-a-model”  turbulence  matrix. 
Coupling  of  individual  effects  will  be  required  to  ensure  the  extent  of  modeling 
capability  as  well  as  the  ability  to  reduce  to  unit  level  problems. 

A  by-product  of  this  Phase  II  effort  will  be  the  development  of  a  “Turbulence 
Tutorial”.  This  tutorial  will  consist  of  a  two-dimensional  boundary  layer  solver  with 
several  turbulence  closure  models  ranging  for  zero-equation  to  two-equation  levels. 
Also  included  in  the  turbulence  turorial  package  will  be  a  database  and  user-manual 
that  will  provide  directions  on  use  of  the  tutorial  to  be  adapted  to  a  PC  that  can  include 
a  post-processor,  graphics,  and  windows  to  allow  students,  researchers,  and  teachers 
the  opportunity  to  learn/explore/teach  turbulence  issues.  While  such  a  product  could 
be  an  invaluable  eductaional  tool  (university),  it  also  provides  a  mechanism  for 
researchers  to  test  turbulence  models  and  variations  of  prior  to  insertion  into  complex 
NS  solvers  (industry  and  government). 
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APPENDIX  A 


MODIFICATION  TO  CEBECI-SMITH  MODEL 

A.1  C-S  Model  with  Compressible  Effects 

For  a  one-dimensional  non-steady  flow,  the  momentum  boundary  layer  equation 
can  be  expressed  as 


du/dt  =  (p/p)  (32n/3y2) 


(A-1 ) 


with  the  boundary  condition  u(o,t)  =  U0  cos  (ft).  The  solution  to  the  above  diffusion 

type  equation  is  easily  obtained.  Cebeci  extended  Stokes  flow  to  include 
compressible  effects  (non-  isothermal)  by  consideration  of  the  non-steady  momentum 
equation  in  the  form 

p  9u/dt  =  dfdy  (|idu/3y) 

(A-2) 

and  introducing  the  transformation  dy  =  pdz  with  the  product  pp“  representing  some 
average  value  over  the  sublayer  region,  Equation  (A-2)  becomes 

du/8t  =  l/(p  p)  (32u/9z2)  (A-3) 


The  solution  to  Equation  (A-3)  subject  to  boundary  conditions  from  the  incompressible 
development,  becomes 

u/uo  =  exp  (-mz)  [cos  ft  -  mz]  / a 


where  m  =  [(f/2)  pp]1/2  =  constant  u+s  (pp/v)1/2  for  u+s  =  (V'P)172.  'rhe  damping 
function  can  then  be  expressed  as 


[1  -  exp  (-yw/  At)) 

where 

y£  =  u£  y  /  vw  and  As  =  26  (y-)1/2(j^)  (A-6) 

Here  one  notes  the  choice  of  the  constant  as  the  incompressible  value  of  26. 
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In  the  above  formulation,  one  could  potentially  use  a  reference  enthalpy  (Eckert, 
for  example)  and  a  viscous  power  iaw  relationship  instead  of  the  average  properties  of 

p  and  p.  In  the  former,  the  equation  of  state  could  be  assumed  such  that  p  ~  1/T 
and  the  latter  as  p.  =  p  (T).  In  this  manner,  direct  use  of  the  Mach  number  and  wall 
temperature  ratio  (Ty^JTa^)  would  be  obtained  to  have  an  explicit  compressible  and 

non-isoenergetic  function.  Most  investigators  have  chosen  to  use  local  values  for  p 
and  p. 


A.2  C-S  Model  with  Pressure  Gradient  and  Mass  Transfer 

Cebeci  recognized  that  the  expression  given  by  Equation  (A -6)  had  to  be 
modified  to  account  for  incompressible  flows  with  pressure  gradients  and  mass 
transfer.  This  was  accomplished  by  redefining  the  parameter  A  in  terms  of  a  sublayer 
friction  velocity  rather  than  its  wall  value,  such  that 

A  =  26  v  „  (xjp  A*  v  JOJp*)1'1 

By  neglecting  the  convective  component  in  the  stream  direction  in  the  momentum 
equation,  the  friction  velocity  in  the  sublayer  can  be  approximated  as 


dzj  dy  -  (vw/v)  x,  =  dp/dx 


(A-7) 


and  for  a  value  of  ys+  =  11.8,  the  solution  of  Equation  (A-7)  allows  for  the  damping 
coefficient  to  become 

A*  =  A*  (vw,  dp/dx) 


The  damping  constant  can  be  expressed  as 
A+  =  26  (pw/p)1/2  (pAw) 


^ compressible  contribution 
-^pressure  gradient  &  mass  transfer  contribution 
incompressible  value 


(A-8) 


where 


N  =  1  -  e*p  ( - 11.8  +  «p  ( - 11.8  (A-9) 

LU  vPw  K  U  tl 
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(A-10) 


P+  = 


Ve  Ue  dl^/d% 
(Xw/pw^3/2 


=  OVpl)  (dp/d%)/(xw/pw)3/2 


—  Vw/lJx 


A.3  C-S  Model  with  Roughness  Effects 

In  the  classic  van  Driest  paper  that  introduced  the  damping  function  for  the  near 
wall  region,  the  concept  was  extended  to  include  rough  wall  effects.  An  artificial  mixer 

(vortex  generator)  was  introduced  at  the  wall  to  mitigate  the  damping  function  (1-e*V/A). 
This  implies  that  the  wall  roughness  encompasses  the  viscous  sublayer  such  that  the 
stabilizing  effect  of  viscosity  is  negated  and  the  effect  of  turbulence  is  felt  all  the  way  to 
the  wall.  Observation  of  experimental  data,  in  law  of  wall  parameters,  shows  the 
velocity  profile  to  decrease  with  increasing  roughness  height,  van  Driest  noted  that 
the  damping  effect  of  the  wall  extends  approximately  y+  =60.  As  a  result,  for  roughness 
elements  which  extend  to  a  value  of  y +  =  60,  the  viscous  influence  of  the  wall  should  be 
nullified. 

In  order  to  accommodate  the  effects  of  roughness  within  the  near  wall  viscous 
region  (i.e.,  kr+  <  60  where  kr+  =  uxkr/vwfor  kr  some  r.m.s  average  height),  a  local 

vortex  generating  factor  is  introduced.  This  is  accomplished  by  adding  a  factor  that  is 
similarand  on  the  same  order  of  the  damping  factor,  namely 

e-yV26  _e  -0  <y726) 

Since  the  vortex  generating  factor  should  be  allowed  to  grow  with  the  size  of 
roughness,  the  factor  was  modified  to  a  form 

exp  (-60  y+f  26  k?-)  (A-11) 

such. that  at  y+  =  kr+  =  60,  the  disturbance  factor  just  offsets  the  damping  factor.  The 
mixing  length  is  then  expressed  as 

tm  =  Ky+  [  1  -  exp  (-yV26)  +  exp  (-60y+/k?)]  (A- 1 2) 


An  examination  of  the  above  shows  that  the  van  Driest  formulation  will  not 
decrease  any  further  once  the  sublayer  has  been  eradicated  by  a  sufficiently  large 

roughness  height  (kr+  >  60)  which  is  in  conflict  with  experimental  evidence 

(McDonald-Fish,  1973).  in  order  to  allow  the  damping  factor  to  be  in  excess  of  unity, 
an  incremental  factor  due  to  roughness  is  introduced,  namely 
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Dkr  =  [  1  +  kr  /  30y+]  exp  (-2.3y+/k?) 


(A-1 3) 


such  that  damping  becomes 


D  =  Ds+ 


(A-1 4) 


Here,  Ds  is  the  original  smooth  wall  damping  term  and  provides  the  sole  means  of 
introducing  the  effect  of  roughness  into  the  prediction  process.  The  mixing  length  is 
now  written  as 


l^ky+D 


(A-1 5) 


Further  comments  concerning  this  development  were  addressed  in  the  text.  It  should  be 
noted  that  the  form  of  Equation  (A-1 5)  is  similar  to  that  which  was  recently  derived  by 
Krogstad  (6-91). 

The  turbulent  viscosity  coefficient  is  obtained  from  the  inner  to  the  outer  layer 
about  a  switch  point  defined  as  follows 


f  Htj@y<yc 
1  Pto  @  y  >  yc 


where  yc  is  the  smallest  value  of  y  where  \iq  =  ty0. 


(A-1 6) 
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APPENDIX  B 


B.1  Baldwin-Lomax  (B-L)  Model 

The  Baldwin-Lomax(1978)  turbulence  model  is  included  here  for  completeness 
inasmuch  as  it  has  become  integrated  into  many  industrial  codes  for  direct  use  or  as  a 
standard  from  which  other  models  are  compared.  Like  the  C-S  model  the  B-L  model 
also  works  well  for  attached  boundary  layer  flows.  However  it  incorporates  features 
that  are  more  advantageous  than  the  C-S  model  by  eliminating  difficulties  in  the  outer- 
layer  length  scale  that  require  a  knowledge  of  the  boundary  layer  edge.  Moreover,  the 
inner-layer  velocity  scale  uses  the  magnitude  of  the  vorticity  vector  as  opposed  to  the 
shearing  strain.  This  provides  3-D  capability  for  flows  where  an  invariant  shearing 
strain  is  not  well  defined.  This  model  has  the  following  features: 


M-d  =  pi2  R  =  p  (k  y  D)2  (cfl{  (B-1) 


where  Ic3|  is  the  magnitude  of  the  local  vorticity  vector  and  the  damping  factor  D  has  the 
same  functional  form  as  in  the  C-S  model.  The  vorticity  term  has  the  same  form  as  C- 
S  when  the  thin-layer  assumption  is  invoked.  Baldwin-Lomax  chose  local  values  for 
the  fluid  properties  in  the  damping  coefficient  of  Equation  (A-8)  and  noted  that  the 
compressible  contributions  become  significant  only  when  large  temperature  gradients 
occur  in  the  near  wall  region. 

The  outer  region  of  the  B-L  model  modifies  the  Clauser  relation  of  C-S  to  give 


jxto  =  0.0268  p  ymax  Fmax  Yo 


(B-2) 


where  y,^  and  correspond  to  the  values  at  the  location  of  the  maximum  value  in 
tha  vorticity  function,  F(y)  =  yDlco|.  The  modified  Klebanoff  intermittency  factor  is 


Yo  -  [  1  +  5-5 


(B-3) 


where  one  notes  that  8  is  replaced  by  ymax/0.3  in  Equation  (5).  The  turbulent 
coefficient  switch  point  is  defined  as  in  Equation  (10).  Table  B.1  provides  further 
details. 
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Table  B.1 .  BALDWIN-LOMAX  TURBULENCE  MODEL  (from  Menter.  6-91) 

The  Baldwin-Lomax  Model 

The  formulation  for  the  eddy  viscosity  is 

/  (dinner  y  -  J /crossover 
\  (ut)  outer  V  >  Vcroasover, 

where  y crossover  is  the  smallest  value  of  y,  for  which  the 
inner  and  the  outer  formulation  are  equal. 

The  inner  eddy-viscosity  is  defined  as  follows: 


/inner 


l  =  ky[\  -ezp(-y+M+)]. 
|w|  is  the  magnitude  of  the  vortidty  and 


The  outer  formulation  is 


(^t) outer  -  KCcpFvta}CeFkleb(v) 


F-wake  =  m*n(ymaxFmax]  C^ymax^i  j I Fmax )• 

(6) 

The  quantities  ymax  and  Fmax  are  determined  fiom  the. 
function 

F(y)  «  yM(l*0  -  exp(-y+M+)).  (7) 

Fmax  is  the  maximum  value  of  F(y)  and  ymax  is  the 
copesponding  y  location.  The  function  Fkiei,{y)  is  the 
Klebanofif  intermittency  function: 

FUet(v)=[l+5j(^i!)6]"1.  (8) 

Hie  quantity  u^*/  is  the  difference  between  the  maximum 
and  the  minimum  total  velocity  in  the  profile.  The  constants 
used  are 

A+  »  2 6&,  Cep  •  1j6;  Cfeleb  =  03; 

Cwk  =  1.0;  k  =  0.4;  K  =  0.0168. 


APPENDIX  C 


C.1  Johnson-King  Model  (J-K):  One-Half  Equation 

The  Johnson  and  King  model  (11-85)  represents  a  nonequilibrium  eddy  viscosity 
model  that  is  formulated  by  an  algebraic  eddy-viscosity  relation  and  an  ordinary 
differential  equation  (ODE)  based  on  a  simplification  of  the  turbulent  kinetic  energy 
field  equation.  The  ODE  is  based  on  the  distribution  of  the  maximum  turbulent  shear- 
stress  which  controls  the  value  of  the  outer  eddy-viscosity  and  ensures  that  lag  effects 
are  accounted  for.  The  model  has  been  successfully  used  to  predict  turbulent 
boundary  layer  flows  including  shock  interactions.  The  model  has  been  updated  by 
Johnson-Coakley  (11-90)  to  include  a  modification  to  the  near  wall  eddy  viscosity 
formulation  as  well  as  incorporating  compressibility  effects.  Menter  (6-91)  has  used 
the  updated  version  of  the  J-K  model  and  compared  results  with  two-equation  models 
as  well  as  the  Baldwin-Lomax  and  Baldwin-Barth  (1-91)  one-equation  model.  Table 
C.1  provides  the  modeling  characteristics  that  include  the  Johnson-Coakley 
modifications  (table  was  structured  from  Menter). 

J-K  attempted  to  develop  a  closure  model  for  a  limited  class  of  flows  (as  opposed 
to  a  universal  model)  which  was  motivated  by  the  ability  of  the  C-S  model  to  perform 
with  known  invalid  assumptions  of  the  turbulent  stress  dependence  on  local  properties 
of  the  flow  (Boussinesq  criteria).  Accordingly,  the  closure  model  without  the 
requirement  to  predict  the  production,  dissipation,  and  diffusion  (as  is  done  for  two- 
equation  models)  was  developed.  For  flows  experiencing  severe  pressure  gradients, 
it  was  assumed  that  characterizing  convection  effects  would  be  essential  with  diffusion 
less  dominant. 

An  ODE  was  developed  that  considers  the  maximum  Reynolds  stress  in  the 
stream  direction  that  prescribes  the  development  of  u‘v‘m  together  with  an  eddy 
viscosity  distribution  across  the  shear  layer  that  is  functionally  dependent  on  uVm 
(subscript  m  refers  to  the  location  of  maximum  trubulent  shear  stress).  It  should  be 
noted  that  this  development  is  in  spirit  to  that  of  Norris-Reynolds  (1975)  and  Wolfshtein 
(1969)  where  damping  of  the  eddy  viscosity  (using  a  turbulence  Reynolds  number 
based  kinetic  energy  characteristic  velocity)  is  not  a  viscous  effect  but  is  the  result  of 

the  near  wall  reduction  of  the  normal  fluctuations  (v')2  due  to  pressure-strain  and  has 
been  verified  by  direct  numerical  simulation  (Rodi,  1-91). 

J-K  provided  the  following  functional  form  for  eddy  viscosity 


Vt  =vto[l  -  exp  (-Vti  /vto)] 


(C-1) 
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where  vti  and  vt0  represent  inner  and  outer  regions  of  the  shear  layer  and,  in  effect, 
can  be  considered  a  two-layer  model.  In  Equation  (C-1),  vt  is  functionally  dependent 
on  vt0  across  most  of  the  shear  layer  (this  is  commensurate  with  production  = 
dissipation  which  is  commonly  used  by  modelers  and  verified  using  DNS  and 
experimental  data).  The  outer  layer  is  used  to  control  the  rate  of  growth  of  -u‘vm.  The 
inner-layer  eddy  viscosity  is  written  as 

Vti  =  k  y  Dj.k  (u/  v'  m)1/2 


where  Dj_k  is  a  damping  term  similar  to  C-S.  The  outer-layer  is  expressed  as 


vto  ■  a  (x)  (0.0168)  Uc8*Yo 

where  a  (x)  represents  an  unknown  modeling  parameter  (it  represents  the  history 
effects  to  account  for  slow  response  of  the  Reynolds  stresses  to  local  changes).  In  the 
above,  the  damping  term  is  given  as 


for 


Dj.k  =  [l-exp(^/Afk)] 


(C-4) 


Ud  =  (u'v'm)l/2  and  Aj.k  —  17 


(C-5) 


It  is  noted  that  the  damping  term  represents  a  modification  of  the  van  Driest 
formulation  to  account  for  y  dependency  instead  of  y2.  The  subscript  m  denotes 
evaluation  when  the  Reynolds  stress  terms-u'v'  reaches  a  maximum.  In  the  outer- 
layer  y0  represents  the  Klebanoff  interm ittency  function. 


Modifications  to  the  J-K  model  were  made  by  Johnson  and  Coakley  (11-90 )  that 
incorporated  a  compressibility  term  to  the  velocity  scale  of  Equation  (C-2)  that  consists 

of  the  density  ratio  (p^p)172.  Length  scales  are  redefined  to  allow  for  the  J-K  model  to 
be  used  for  flows  that  are  less  severe  than  those  experienced  in  separation.  Details  of 
these  modifications  are  found  in  Table  C.1  and  a  discussion  of  same  in  the  Johnson- 
Coakley  paper. 
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Table  C.l.  JOHNSON-KING  TURBULENCE  MODEL  WITH  JOHNSON-COAKLEY 

MODIFICATIONS  (from  Menter,  6-91) 


The  Johnson-King  Model 

The  JK-modd  is  a  nonegnilibrinm  eddy-viscosity 
model,  which  is  composed  of  an  algebraic  eddy-viscosity 
formulation  and  an  onhnary  differential  equation  (ODE). 
The  ODE  is  a  simplified  Reynoids-stress  equation  along 
the  path  of  the  maximum  turbulent  shear-stress.  It  enw|fw*g 
the  value  of  the  outer  eddy-viscosity  i /to  and  ensures  that 
history  effects  are  property  accounted  far. 

The  model  has  recently  been  modified  to  account 
for  a  thin  logarithmic  region  close  to  the  wall  which  exists 
even  in  the  presence  of  strong  adverse  pressure  gradients. 
This  new  formulation  was  used  for  the  prerent  investiga¬ 
tion  (note  that  the  compressible  formulation  depends  on  the 
variable  density). 

The  eddy  viscosity  vt  has  the  following  farm: 


i/t  =  i/totanh(i/ti/i/to) 


with  vti  and  vto  defined  as  follows: 

vti  =  zAeyuj 
vta  =<r(x)  (0.0168ue**7) 

In  this  formulation  6*  is  the  displacement  thickness  of 
the  boundary  layer,  tte  is  the  velocity  at  the  edge  of  the 
boundary  layer  and  7  is  Klebanaff’s  interim  tiency  function: 
7  =  • The  vclocity  scale  n,  which  appears  in 

the  formulation  of  vti  is  a  Mending  between  the  friction 
velocity  xir  and  um  =  (— ttVm)1/2^  index  m  refers 
to  the  location  of  the  maxrnmm  turbulent  shear  stress): 

u«=ttr(l-72)  +  UmTl 

72  =  tanh(y/Z/c) 


L'  = 


«r 


- Vm 

u r  +  Um 

The  damping  function  D  is  defined  as  follows: 


D- 


where  A+  =  17  and  up  =  moifttm,  ur). 

The  function  o(x)  is  determined  from  the  solution  of 
the  ODE  for  the  maximum  Reynolds  shear  stress: 

(-Wm)'/2  = 

oi(-u't ;'m)  ax 

- is—  On 

(— tt'l/m) 

In  this  equation  Lm  is  the  dissipation  length  scale,  aj  the 
ratio  of  Reynolds  shear  stress  to  turbulent  kinetic  energy, 
Dm  the  turbulent  diffusion  and  (-uVm)e?  the  maximum 
Reynolds  shear  stress  when  convection  and  diffusion  are 
zero.  (— tiVm)e<7  can  be  obtained  by  setting  a(x)  =  1 
in  Eq.  (  )  and  replacing  %  in  Eqs.  (  )  and  (  )  by 
(u m)eq-  Lm  and  Dm  are  defined  as  foUows: 


Lm  *  0-4ym*  Vm/ 8  £  0.225 
Lm  =  0.095,  y m/8  >  0225 


Dm 


Cdif(-U'v'm)2/2 
aj6[0.7  —  (y/5)m) 


|1  -  <r{x)lf2\ 


The  value  of  o{x)  is  chosen  for  every  profile  in  a  way,  that 
ensures  that  the  maximum  shear  stress: 


-(uVm)oct  =  h(^  +  ^)]m 


is  equal  to  the  maximum  shear  stress  (— uVm)o D E  com* 
puted  from  the  ODE,  Eq.  (  ).  This  is  achieved  iteratively 

by  _ 

a(*r+'=<Kx)w(7“^l00E 

(-trt fm)aet 

Only  one  iteration  is  performed  per  time  step.  The  con¬ 
stants  are: 
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O!  =025 
Cdif  -  05 
Cdif  =  0.0 


<r(x)  >  1 
a(x)  <  1 


APPENDIX  D 


D.1  Norris-Reynolds  Model  (N-R) 

The  authors  chose  the  turbulent  kinetic  energy  closure  concept  of  Prandtl  for  the 
eddy  viscosity  such  that 

vt  =  vt  (k,  t) 

where  k  is  prescribed  by  the  turbulent  kinetic  energy  transport  equation.  The  length 
scale  (i)  is  algebraically  specified.  In  addition,  a  closure  relationship  for  the 
dissipation  rate  is  required  such  that  e  =  e  (k,  i,  v). 

The  dissipation  rate  e  which  is  a  length  scale  in  the  eddy  viscosity  relation  of 
Equation  (11)  and  source  a  term  in  the  kinetic  energy  field  equation,  is  determined 
from  the  relation 


where 


vt  =  k112  v 


(D-1) 


V  =  Q  y  [1  -  exp  (-C3  Rt)] 


(D-2) 


for 


Rt  =  P  km  y/p 


(D-3) 


The  eddy  viscosity  that  is  proportional  to  a  turbulence  Reynolds  number  is  that 
proposed  by  Prandtl  and  asserts  that  the  same  large  eddies  responsible  for  turbulent 
mixing  are  also  those  containing  the  bulk  of  the  turbulent  motion.  The  constants  C^,  Ct 

and  C3  are  to  be  determined.  The  one  major  difference  in  the  damping  term  of 
Equation  (D-2)  is  the  use  of  kinetic  energy  velocity  as  opposed  to  the  wall  shear 
velocity  in  the  van  Driest  relation.  It  is  believed  that  the  use  of  the  kinetic  energy  as  a 
velocity  scale  is  one  criteria  for  the  success  of  simple  turbulence  models  to  predict 

severe  pressure  gradient  flows  (since  it  has  been  demonstrated  that  ux  is  no  longer 
linear  for  these  conditions). 

Norris-Reynolds  and  Wolfshtein  (1968)  allowed  for  the  pressure  velocity  term 
and  the  triple  velocity  correlation  term  to  be  modeled  in  a  gradient  diffusion  manner 
such  that 
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—  ( k  u'i  u'i  u'j  +  p' u'i  Sij/p')  =  —  (v  dk/3xj) 
dxjV  '  3*j 

For  this  condition,  the  source  term  can  be  expressed  as 


e  =  CDk3/2Ae 


(D-4) 


which  is  based  on  dimensional  analysis  and  is  a  result  found  by  Kolmogorov,  i.e., 
e  oc  k372/!.  The  outer  layer  (of  the  boundary  layer)  should  be  controlled  by  the  factor 
CD  V?,2K  which  will  require  modification  to  mitigate  the  size  of  the  eddies  by  allowing 
for  the  viscosity  to  be  a  dissipative  agent  for  the  largest  eddies. 

Variations  on  the  length  scale  have  been  introduced.  Woifshtein  recommended 
a  form  similar  to  Equation  (D-2),  namely 

h  =  Q  y  [  1  -  exp  (-Rt/Ae)]  (D-5) 

while  Norris-Reynolds  chose 

ie  =  Qy(l  +  C4/RT)4  (D-6) 

and  meets  the  requirement  that  dissipation  at  the  wall  is  non-zero.  The  dissipation 
becomes  the  Kolmogorov  function  in  regions  away  from  the  wall  and  has  a  functional 

relation  e  <*=  vk/i2  in  the  near  wall  region.  Here,  one  notes  that  i  «=  y  near  the  wall  as 
opposed  to  the  van  Driest  equation  which  implies  i «  y2. 

The  purpose  of  the  van  Driest  damping  factor  is  the  near  wall  suppression  of  i  in 
order  to  make  the  eddy  viscosity  vt  =  i2  I9u/9yi  very  small  in  the  sublayer.  Accordingly, 
Equation  (D-1)  was  developed  in  place  of  the  van  Driest  damping  factor  for  near  wall 
dissipation.  Here  the  exponential  near  wall  damping  factor  was  removed  from  the 
length  scale  requirement  and  became  a  direct  part  of  the  eddy  viscosity  specification. 
Moreover,  the  small  eddy  characteristics  of  the  near  wall  region  are  retained  while  the 
proper  sublayer  modeling  of  the  dissipation  was  achieved  and  provides  a  mean  flow 
model  that  is  valid  all  the  way  to  the  wall. 

In  the  above,  two  distinct  lengths  have  been  identified.  A  length  scale 
appropriate  to  viscosity  (O  and  a  length  scale  characteristic  for  dissipation  (ie).  The 
constants  used  in  the  length  scales  are  as  follows: 

is  same  as  k- e  (  =  0.09) 
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Cv  is  chosen  for  conformity  to  the  log  law 


and 


(D-7) 

C3  =  26/Ap.  A+ 

(D-8) 

where  A+  is  equal  to  26  (van  Driest)  and  50.<  An  <  70.  Again  A+  can  include  terms  to 
account  for  surface  roughness  and  mass  transfer  as  well  as  pressure  gradient 
conditions.  The  remaining  constants  have  been  selected  as: 


CD=  unity 

C4  =  5.3  (D-9) 

Rodi  (1-91)  noted  that  damping  due  to  the  exponent  of  Equation  (D-2),  as 
described  above,  is  not  a  viscous  effect  but  is  a  result  of  the  near  wall  reduction  due  to 
pressure-strain  that  has  been  confirmed  by  direct  numerical  simulation  (DNS). 
Moreover,  damping  of  the  dissipation  length  te,  Equations  (D-5)  and  (D-6),  are  not  in 
agreement  with  DNS  data  which  suggests  that  ie  =  C^y;  however,  the  multiplication 
by  the  damping  term  provides  good  engineering  estimates  for  wall  conditions. 


APPENDIX  E 


E.1  Goldberg  Model  (6-91) 

Recognizing  that  large-scale  eddies  are  characterized  by  a  development  rate 
much  different  from  small  scale  dissipative  rate,  a  one-equation  model  was  developed 
using  two-time  scales  of  turbulence.  These  time-scales  provide  length  scales  for 

diffusion  (i^)  and  dissipation  (ie)  that  are  similar  to  those  developed  by  Norris  and 

Reynolds  or  Wolfshtein.  The  velocity  scale  is  determined  by  solutions  to  the 
turbulence  kinetic  energy  PDE  (standard  form  of  k-e).  The  time  scales  for  the  large  and 
small  eddies  (using  the  Kolmogorov  scale)  are  given  by 


tfc  -  k/e ,  Large  Scale 
tk  -  Vv/e,  Small  Scale 

where  x>  is  the  kinematic  molecular  viscosity. 


(E-1) 


The  dissipation  is  algebraically  expressed  in  terms  of  k  and  a  length  scale 
throughout  the  boundary  layer  as  provided  by  Wolfshtein  i.e.  Equations  (D-4)  and 
(D-5)  with  a  condition  on  the  y  coordinate;  namely 


y=  min  { y,  Cjy5/ic}  (E-2) 

where  5  is  the  boundary  layer  thickness  and  Kthe  von  Karman  constant. 

Combining  Eqs.  (E-1)  with  the  velocity  scale  k,  the  eddy  viscosity  associated 
with  the  time  scales  are 

vk  ~  ktic=  Cjx  k2/e ,  see  eq.  (I  \) 


and 


ve  ~  k  te  =  Ce  k  Vv7e  [1  -  exp  (-A^Rt)]  (E-4) 

where  a  damping  term  has  been  added  to  Eq.  (E-4).  The  constant  C£  is  obtained  from 
comparison  to  the  N-R  results  i.e.  Eqs.  (D-1)  and  (D-2)  such  that 

Q  =  VIkCJ/4  (E-5) 

Goldstein  combined  Eqs.  (E-3)  and  (E-4)  into  a  single  expression  which  has  the 
correct  behavior  in  the  near-wall  region  as  well  as  away  from  the  wall,  such  that 
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Cn  fn  =  (C^  +  Q  VvE/k)  [1  -  exp  (-A^  Rt)]  (E-6) 

where  has  limits  pending  on  Rj»1  and  Rj  «  1.  This  led  to  an  eddy  viscosity 

expression  for  the  entire  layer  given  as  vE  =  k2/e  which  is  the  well  accepted 

constitutive  equation  for  k-  e  models  using  wall  damping  functions.  What  makes  the 
Goldstein  work  different  from  the  other  one-equation  models  is  the  definition  of  the 
length  scales.  Using  the  algebraic  definition  for  the  length  scale  (Eqs.  D-4  and  D-5) 
together  with  the  velocity  scale  (k1/2)  and  Eqs.  (E-3)  and  (E-4),  the  length  scales 
corresponding  to  large  and  small  scale  eddies  become 

ik  =  *  qi/4  y  [1  -  exp  (-Af*  RT)]  (E-7) 


le  =  [^*^1/2  y  (1  -  exp  (-Ae  Rt)]  m  [1  -exp  (-A^  RT)]  (E-8) 

where  provisions  on  y  are  specified  by  Eq.  (E-2).  Arguments  are  presented  relative  to 
scaling  criteria  of  Eq.  (E-8)  as  y  — >  o.  In  particular,  the  behavior  of  the  turbulence  KE 
in  the  viscous  sublayer  that  provides  a  density  ratio  with  the  characteristic  velocity 
scale.  This  was  also  shown  by  Gorski  (1-86)  as  well  as  by  the  Johnson  and  Coakley 
(11-90)  correction  to  the  velocity  scale  of  the  Johnson  and  King  model.  Because  of  the 
introduction  of  a  simple  multi-scale  procedure  and  its  similarity  to  models  that  have 
had  demonstrated  success,  this  model  should  be  considered  for  further  evaluation. 
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APPENDIX  F 


F.1  Baldwin-Barth  Model  (B-B)  :  One  and  One-Half  Equation  Model 

The  B-B  turbulence  model  (1-91)  considered  a  transformation  of  the  basic 
physical  variables  in  the  k-e  transport  equations  in  order  to  eliminate  the  numerical 
difficulties  in  the  near  wail  region  associated  with  these  two-equation  models.  Table 
F.1  provides  details  of  the  turbulence  structure  as  summarized  in  the  Appendix  of  the 
B-B  paper.  The  authors  chose  a  field  equation  for  a  turbulence  Reynolds  number 
based  on  an  eddy  viscosity  concept  for  high  Reynolds  number  flow  with  the 
distribution  away  from  the  wall  governed  by  Cji  ^2/e  [combines  Equations  (D-1) 
and  (D-4)];  in  particular,  the  turbulence  Reynolds  number  Rey  is 


Rej  =  k2/oe  (F-1) 

The  above  is  substituted  into  the  standard  k-  e  equations  and  the  results  are 
combined  to  provide  a  one-equation  model  that  encompasses  turbulent  kinetic  energy 
and  its  dissipation.  Equation  (F-1)  is  rearranged  to  provide  two  regions  in  the  shear 
layer  consisting  of  an  outer-layer  (k2)where  production  is  equal  to  dissipation  and  an 
inner-  layer  (k^where  production  requires  damping;  namely 

e  =  k2/u  Rej  =  (k1  +  k2)Aj  Rej  (F-2) 

For  the  outer  region  p  ■  e,  and  k1  »  k2,  then  k2  =  x>  Reyp,  and  the  turbulent  kinetic 
energy  becomes 

k=  (uRey  p)1/2  (1  +  kg/ty  (F-3) 

When  Equation  (F-3)  is  substituted  into  the  standard  high  Reynolds  number  k-e 
equations,  there  results 

D  (v  Rct)  /Dt  =  (C2-  Ci)  (v  Rex p )m+  ( v  +  )  V2(v  Rer )  - 

\  <Jk  /  (F-4) 

^(Vvt)*V(vReT)-(2-C1)j^_(vRer  P)-(2  -C2)k2 

The  system  is  closed  by  substituting  Equation  (F-3)  into  the  kinetic  energy  equation  of 
the  k-e  equations  to  give 

Dk/Dt  =  -2  (p/vRer)1'2  k2  -  k^Rer  +  V»(v  +  ^)Vk  (F-5) 
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At  this  junction,  the  B-B  approach  follows  a  path  similar  to  thehybridr  modelers 
(Chen  and  Patel,  11-88),  (Rodi,  1-91).  Here,  Equation  (F-4)  is  solved  for  k2  ~>  0  for 
solution  in  the  high  Reynolds  number  region  (which  is  most  of  the  shear  layer)  and 
develops  a  length  scale  to  accommodate  damping  in  the  near  wall  region.  One 

should  keep  in  mind  that  since  field  equations  for  k-e  were  combined,  the  results 

represent  solutions  to  the  k-  e  two-equation  model  for  the  regions  away  from  the  wall. 
To  accommodate  the  near  wall  region,  the  turbulence  Reynolds  number  is  split  into 
two  factors 

Rd  =  Rer  h  (Rei) 


where  f3  is  a  damping  function  that  allows  Rej  «  Rej  for  large  Reynolds  numbers. 
Using  the  damping  functions  similar  to  those  of  van  Driest  or  N-R,  there  results 

vt  =  vC,f,Rer  =  vC,f,f3Rer  (F-7) 


and 


k2/v  Rer  =  (ki  +  k2)2  /  v  f3  Rer  (F-8) 


where  k.,2  =  u^Rej.  Using  the  thin-layer  assumption  where  p  -  ut  (Uy)2,  the  field 
equation  becomes 

DReT/Dt  =  (C2f2-C1)(C^fHf3)1/2Rer  — +  (v  - 

By  dy2  (F-9) 

(l/crE)  (9v t/9y)  B  (Rer)  /  9y 


which  becomes  the  basis  for  the  one-equation  model  defined  in  Table  F.l.  The 
damping  term  f3  has  to  be  selected  to  allow  for  k1  to  be  the  dominate  part  of  k  in  the 

near  wall  region.  Details  for  the  selection  of  the  damping  function  are  provided  in  the 
B-B  paper. 


Menter  (6-91)  compared  the  B-B  model  to  adverse  pressure  gradient 
experimental  data  as  well  as  other  turbulence  models.  It  is  not  clear  what  velocity 
scale  was  used  in  the  damping  terms.  It  appears  that  the  van  Driest  shear  velocity  was 

used,  i.e.,  y+  instead  of  Rej  where  Ux  is  used  in  y+  and  k1/2  is  used  in  Rej.  It  was 

found  that  displacements  (hence,  pressure)  were  predicted  with  good  accuracy,  but 
could  not  predict  profile  shapes  in  either  the  inner  or  outer  shear  layer  regions. 
Moreover,  wall  shear-stress  distributions  are  similar  to  ones  produced  by  the  B-L 
model. 
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Table  F.l.  BALDWIN-BARTH  ONE-EQUATION  TURBULENCE  MODEL 

(from  Baidwin-Barth,  1-91) 


=  (  _  )Sl  +  (l/  +  ±l)Gi v\urt) 

Dt  <TR 

In  this  equation,  we  use  the  following  functions: 


—  — (c«3  —  c.Jv^T/zc2 
<r« 

<tr  =0.25  <r« 
ut  =c,1{uRt)D\D2 

fit  =put 

Di  =1  -  exp  (-y^/A+) 
D2  =1  -  exp  (-y+/Af) 


P=ut 


dUj  dU: 
dxj  ^  dx 


±\dVi_2i/  fdUk 
i  )  dxj  3  ‘  \  dxk 


My* )  =r1  +  (i-^-)(-^  +  d,d2)(  ■Jd^ 

c«j  c«a  KV 

eip(-y+M+)  D* + 


Gi  =1  + 


(F*+*o/B7) 

(w4  +  B^Fjt  -f-  «o) 


„  (VvRr)  •  (VvRt)  k 2,  , 

- - - - — M 

i/itr  c»* 


For  all  calculations  we  have  used  the  following  constants: 


k  =0.41,  c«t  =  1.2,  c«3  =  2.0 

Cp  =0.09,  A+  =  26,  At  =  10 
Bi  =0.4,  B2  =  0.01,  co  =  10“10 


We  also  recommend  the  following  boundary  conditions  for 

1.  Solid  Walls:  Specify  Rt  —  0. 

2.  Inflow  (V  •  n  <  0):  Specify  Rt  =  (Rt)oo  <  1. 

3.  Outflow  (V  •  n  >  0):  Extrapolate  Rt  from  interior  values. 


— (V*,)  ■  V(vfiT) 

<7t 


*Mt)  Dt)) 
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Since  the  B-B  model  is  relatively  new  and  untested,  further  work  will  be  required 
prior  to  a  full  assessment.  The  N-R  damping  functions  should  be  considered  for  the  B- 

B  model  which  incorporate  kinetic  energy  to  the  1/2  power  (k1/2)  and  uses  a  modified 
Kolmogorov  function  (for  length  scale)  relating  kinetic  energy  dissipation  that  is 
applicable  throughout  the  shear  layer. 
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q.  =  -  [u/Ptl  -  uT/PrT] 


ah 

3x, 


(  ") 


Wilcox  k -to  Turbulence  Model 
(from  Wilcox,  6-91) 


For  general  compressible  turbulent  fluid  flow,  the  complete 
sets  of  eauauons  that  constitute  the  Wilcox  k-u  two-equa- 
uon  model  and  muitiscaie  model  are  written  in  terms  of 
Favre  mass-averaged  quantities  as  follows. 


Conservation  of  Mass,  Momentum  and  Energy 


Bo  3 

— - —ion.)  =  0 

3t  ox. 


3  3  3  r  c  -  , 

—  (pu,)  -  —ipu.u,)  =  r-f-pS,,  -  r.J 
0 1  ax.  ox. 


(  1) 


(  2) 


3  3  3  r  3k  . 

—  <pE)*  —  <Pu£)  =  f^-u,  -  q,  +  iu~o"uT)—] 

A#  Aw  *  Av  1  1  r*  V 


3*, 


3x, 


ox. 

t 

( 


3) 


where  t  is  time,  x,  is  the  position  vector,  u.  is  velocity,  p  is 
density,  p  is  pressure,  u  is  molecular  viscosity.  ?.,  is  the  sum 
of  molecular  and  Reynolds  stress  tensors,  and  a  is  the  sum 
of  molecular  and  turbulence  heat  flux  vectors,  fn  Equation 
(  3).  the  quantities  E»e^k+Viu,u.  and  H=h«-k  +  y;u,u,  are 
total  energy  and  total  enthalpy,  respectively,  with  h =e  +  p/p; 
e  and  h  denote  internal  energy  ana  enthalpy.  Also  5,,  is  the 
Kronecker  delta,  and  k  is  the  turbulence  kinetic  energy  that 
is  determined  by  the  following  equations. 

Turbulence  Kinetic  Energy,  Specific  Dissipation  Rate 
3  3 

—  (pk)  +  — (PUjk)  *  r  — - Bmou k 

3t  3x,  r  3x. 


3  .  3k, 

3  3  « 


(  4) 


-  apui(2QmaQwm)* 


B  .  3u , 


5) 


The  quantity  <a  is  specific  dissipation  rate.  r„  is  Reynolds 
stress  tensor,  aqd  pT  is  eddy  viscosity.  The  parameters  a.  3. 
0 •.  a ,  a*  and  (  are  closure  coefficients  whose  values  are 
given  below.  We  dose  the  system  of  equations  as  follows. 


Constitutive  Relations 
1 

h  -  2q(Sij  -  - +  r„ 


(  6) 


ur  =  a*pk/ia 


1  3u,  3u. 


CO 


=  €/K 


1  3u,  3u. 


ox, 


(  8) 


(  9) 


The  quantities  Pr,  and  PrT  are  laminar  and  turbulent 
Prandtl  numbers,  respectively.  and  aa  is  an  additional 
closure  coefficient  that  is  specified  below.  We  need  just 
one  additional  constitutive  relation  for  the  k-u  model. 
Specifically,  we  assume  the  Reynolds  stress  tensor  is  pro¬ 
portional  to  the  mean  strain-rate  tensor.  The  resulting 
relationship  is  as  follows. 

T wo- Equation- Model  Reynolds  Stress 


1  3uk  2 

r.  =  -"rfS.,  -  j  — «„]  •  JOW,, 


(  10) 


The  multiscale  model  computes  each  component  of  the 
Reynolds  stress  tensor  separately.  The  model  introduces 
two  enerey  scales  corresponding  to  upper  and  lower  parti¬ 
tions  of  tne  turbulence  energy  spectrum.  The  quantity  pT 
denotes  the  upper  partition  contribution  to  the  Reynolds- 
stress  tensor  whiie  ee  denotes  the  energy  of  the  eddies  in 
the  lower  partition.  The  additional  equations  needed  for 
closure  are  as  follows. 


Multiscale- Model  Reynolds  Stress 


r,i  ^  pe3:i 


a7<pT:>)  ‘  a7,WkT")  =  *P-1*E-I 


8\ 


(  11) 


(  12) 


3  3  .  »  3u, 

—  (p(k-e)i  «■  —  (pu.lk-e)]  =  (l-a-fllr,,— 


3t 


ox, 


3.x, 


-  3*Ptak(  1  -  e/'k)3  : 

E,,  =  -Ctd*cj{rt|  *  r  Pk5tJ]  *  aP,,  *  dD*, 

*  +  r<5‘«ukn-e/k)3 


(  13) 


I  14) 


P,,  =  7,r 


3Uj  3W  3w„  3u™ 

3xm  *  **  3x»  '  D,j  *  3x,  *  Tj"  3x,  (  15) 


Finally,  these  equations  involve  a  number  of  closure  coeffi¬ 
cients  whose  values  are  given  by  the  following  equations. 


Closure  Coefficients 
B  =  3/40.  Bm  =  9/100.  o  *  1/Z  a*  =  1/2 
am  -  1.  a  ■  [tf(l  +  iVJ1r)-OK2vT^]/tf• 
a- 42/55.  B  -6/55..  9  -  1/4  ,  |  -  l 


(  16) 
(  17) 
(  18) 
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k-e  Turbulence  Model 
(from  Speziale  et  al.f  6-91) 


For  simptierty,  ire  will  lauict  our  analysis  to  i 
prrasihiw  turbulent  flows  (howtvtr,  the  crucial  conclu¬ 
sions  that  will  be  drawn  carry  over  to  fmmprrssihlw 
flows).  The  mean  velocity  TX  and  mean  prescnre  f  are 
solutions  of  the  Reynolds  averaged  Navier-S  takes  and 
continnity  equations  given  by 


at  a*  a*i 


an* 

Bxt 


=  0 


(1) 

(2) 


where  Tjy  =  — ujuj  is  the  Reynolds  stress  tensor,  v  is 
the  kinematic  viscosity,  and  the  usual  Einstein  sum¬ 
mation  convention  applies  to  repeated  indices.  We 
will  consider  the  commonly  used  two-equation  mod¬ 
els  based  on  an  eddy  viscosity  where 


*»/  =  +  **r 


(3) 


and 


(4) 


given  that  K  s  is  the  turbulent  kinetic  energy, 
e  =  vdv^/ dxj  is  the  turbulent  dissipation  rate, 
and  Cft  is  a  dimensionless  constant  at  high  turbulence 
Reynolds  numbers.  In  two  equation  models,  transport 
equations  are  solved  for  any  two  linearly  independent 
variables  constructed  Rom  K  and  c.  In  the  K  —  e 
model,  modeled  transport  equations  fin  K  and  c  are 
solved;  in  the  K—u  model,  modeled  transport  equa¬ 
tions  fin  K  and  the  reciprocal  turbulent  time  scale 
u  =  c/K  are  solved;  and  in  the  K —r  model,  modeled 
transport  equations  fin  K  and  the  turbulent  time  scale 
rsK/e  are  solved.  The  exact  transport  equations  fin 
K  and  e  areas  follows  [8]: 


JDJL  bin  _  _}  _ 

-  ■Sr-’**—- V+“VK 

where  D/Dt  =  d/dt  +TT-V.  In  (5)  -  («), 

V,  =  to-L  lit**} 

axtyaxjdx,) 

a*f  dmh  a»h  J 


(8) 

(«) 

(7) 

(«) 


are  turbulent  transport  terms,  and 


V. 


_  2,55^ 

day  a%j  fist  fis<  dsk  da* 

,2,555_2 

dxk  a*j  dxj  Bzj  dzjdxi 


t») 


B  S 

axfdxk  dxjdxk 


(10) 


art,  respectively,  the  production  and  destruction  of 

rli««iyittUn  ^fpwo 


It  is  a  straightforward  matter  to  show  that  near  a  wall 
K  =  0(y*),  c  =  0(1),  r  =  0(y*)  (14) 

^  *  0(1),  =  ©(y3),  *7  =  ©(y3)  (15) 

**  =  0(y4),  **  =  ©(y3),  V  =  ©(y3)  (16) 

1>  =  O(y),  7>.  =  O(y),  *.  =  0(1)  (17) 

D,  =  0(1),  VaJT  =  0(1),  V2e  =  0(1)  (18) 

where  V  S  rj/dn</dxy  is  the  turbulence  production. 

An  asymptotic  analysis  of  the  JT  —  e  model  will  be 
conducted  first.  In  the  K—t  model,  the  eddy  viscosity 
near  a  wall  is  taken  to  be  of  the  form 

**  =  •  (19) 

The  asymptotic  analysis  presented  in  this  — in¬ 
dicates  that  /„  =  0(l/y)  near  the  wall  since,  dne  to 
(15),  vr  mart  be  of  ©(y3)  in  tins  region.  Of  course, 
sufficiently  far  from  the  wall  /„  assumes  a  value  of  1. 
( Cp  is  a  constant  which  is  typically  taken  to  be  0.09). 
The  turbulent  transport  term  P  in  the  kinetic  energy 
equation  (5)  is  modeled  using  a  gradient  transport  hy¬ 
pothesis: 
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Situ  and  Schetz,  (6-91)  considered  compressibility  effects  from  the  Reynolds 
stresses  for  thin  shear  layers.  Here  shear-stress  is  given  as 


x'xy  =  -p  u'  v'  -  u  p'  v'  -v  p'  u'  -p'  u'  v' 

allowing  for  u  >  >  u  and  the  neglecting  the  triple  correlation,  the  shear  becomes 

T'xy  =  -p  uV-up'v'  0*2) 

Shang  (7-74)  considered  the  first  and  last  terms  on  the  RHS  of  eq.  (1-1)  to 
introduce  compressibility.  The  impact  of  the  triple  correlation,  as  reported  by  Marvin 
Coakley  (11-89),  had  no  effect  up  to  Mach  5,  but  showed  some  effect  at  higher  Mach 
number  conditions. 


Situ  and  Schetz  provided  for  a  length  scale  for  the  scale  for  the  second  term  on 
the  RHS  of  equation  (1-2)  which  together  with  gradient  transport  theory  establishes  an 
effective  turbulent  Schmidt  number  for  mixtures  or  a  turbulent  Prandtl  number  for 
homogeneous  gases.  The  Reynolds  stress  terms  in  equation  (1-2)  are  given  as 


for 


-  p  u'  V  -  p  l£,  (du/dy)2 

(1-3) 

—  2  u  dp  du 

-  u  pv  =  im  - - 

d-4) 

s  dydy 

S  ~  Cu  iu/Cp  ip 

(1-5) 

=  Constant  •  Cu  i2 

In  the  above,  one  recognizes  the  usual  definition  of  the  Prandtl  concept. 
Equation  (1-3),  and  the  introduction  of  the  new  length  ip  due  to  compressibility.  For  thin 


layers 


Ht  =  lm|P 


9u+1Lap 

dy  s  By 


-p'v  =Dt 


where  from  gradient  transport  theory  it  is  shown  that 

' _  d  =Jk=l2  JL^H. 

1 3y  Set  m  Set  ay 

dp  du 


and 


-upV=£  £ 


SCt  ty  &y 


0-6) 


d-7) 


(1-8) 
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If  the  turbulent  Lewis  number  is  unity,  S  =  Prt  and  one  can  allow  for  S  =  Prt  =  Scr 

=  0.8. 

The  authors  compared  their  work  to  six  test  cases  that  include  compressible  flow 
over  a  flat  plate  (M  =  4.67  and  6.55),  tangential  slot  injection  (sonic  jet  and  supersonic 
jet)  and  shock/boundary  layer  interacting  flows.  While  the  new  model  does  not  invoke 
ad-hoc  techniques  and  success  was  shown  to  a  range  of  test  conditions,  it  is  not  clear 
what  form  of  mixing  length  characteristics  were  used  by  the  authors.  Recognizing  use 
of  the  original  Prandtl  concept  as  stated  by  the  authors,  they  did  not  show  if  damping  was 
used,  if  the  length  scale  was  proportional  to  y  of  y2,  or  what  form  of  characteristic  velocity 
scale  was  used,  i.e.,  shear  velocity  or  a  form  of  the  kinetic  energy.  An  inspection  of  the 
original  AIAA  paper  (89-1821)  did  not  reveal  any  further  insight. 

Since  the  method  is  new  and  has  not  been  tested  in  the  literature  by 
independent  sources,  this  work  will  require  further  evaluation.  The  use  of  the  double 

length  scale  which  is  also  provided  by  the  hybrid  modelers  (i.e.,  y  and  i  for  viscous 

and  diffusion  effects)  and  the  ability  of  the  model  to  handle  the  wall  jet  problem  warrant 
further  investigation. 
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